Intrusion and extrusion of liquids in highly confining media: bridging fundamental research to applications

ABSTRACT Wetting and drying of pores or cavities, made by walls that attract or repel the liquid, is a ubiquitous process in nature and has many technological applications including, for example, liquid separation, chromatography, energy damping, conversion, and storage. Understanding under which conditions intrusion/extrusion takes place and how to control/tune them by chemical or physical means are currently among the main questions in the field. Historically, the theory to model intrusion/extrusion was based on the mechanics of fluids. However, the discovery of the existence of metastable states, where systems are kinetically trapped in the intruded or extruded configuration, fostered the research based on modern statistical mechanics concepts and more accurate models of the liquid, vapor, and gas phases beyond the simplest sharp interface representation. In parallel, inspired by the growing number of technological applications of intrusion/extrusion, experimental research blossomed considering systems with complex chemistry and pore topology, possessing flexible frameworks, and presenting unusual properties, such as negative volumetric compressibility. In this article, we review recent theoretical and experimental progresses, presenting it in the context of unifying framework. We illustrate also emerging technological applications of intrusion/extrusion and discuss challenges ahead. Graphical Abstract


Introduction
Intrusion/extrusion of liquids in porous systems, i.e. the process of a liquid penetrating into a previously vapor/gas-filled cavity bounded by some walls or the opposite process of the liquid exiting from the confining environment letting vapor/gas to occupy the vacuum left behind, is important for many technological applications, among which are the separation of liquids [1,2], liquid-phase chromatography [3,4], energy damping, conversion, and storage [5,6], porosimetry [7][8][9], biological and bioinspired channels [10][11][12][13], and many more. Indeed, intrusion/extrusion can also be used to investigate fundamental and elusive aspects of multiphase solid/liquid/phase systems, such as the dependence of the surface tension on the curvature of the meniscus [14][15][16][17], i.e. how much energy it costs to create a liquid surface of prescribed area beyond the simple planar case, or the line tension, i.e. the energy cost of putting the three phases in contact along a line of unit length [18].
Before we continue our discussion, it is useful to remind ourselves that different communities employ different terminologies to refer to similar phenomena or sometimes use similar terminology to refer to different phenomena. In the paragraph above, and in the following, we frequently use the terms intrusion and extrusion to refer to a liquid entering and exiting from a confining environment, a pore, a cavity, regular, or irregular textures of solid surfaces. In some communities, intrusion and extrusion is synonymous of wetting and drying pores and cavities. However, the equivalence between these terms is neither obvious nor universally accepted. For example, in the physics community, the term wetting [19][20][21] is used to describe a surface phase transition in which a macroscopic amount of one phase, typically the liquid, is adsorbed at a single, attractive wall, while a different phase, the gas, is stable in the bulk and sufficiently far from the wall (at ambient conditions). In this line of reasoning, drying can be considered as wetting of a single repulsive wall by the gas phase. However, under suitable external conditions, e.g. tensile condition or overpressure, attractive walls can be dried and repulsive ones can be wetted. In other communities, wetting and drying of a single wall might be used to refer to the local increase and decrease of density compared to the bulk value at an attractive and repulsive wall, respectively. Remarkably, this association between wetting/drying and attractive/repulsive walls must not be confused with hydrophobicity/hydrophilic, which is connected to density fluctuations of the liquid rather than the local absolute value of the quantity.
It might also be useful to know that the wetting or drying of a pore, intrusion, and extrusion as discussed in this review is also referred to as capillary condensation or evaporation [22]. While for a long time, there was a clear distinction between wetting and capillary condensation or between drying and capillary evaporation in the physics literature, more recently it was shown that by a smooth geometrical transformation of the planar wall into a (capped) capillary, the distinction between these phenomena is more subtle than expected [23].
In this review, we will also use the terms 'hydrophobic' and 'hydrophilic'. In some communities, these terms are used loosely to identify a wall or material repulsive or attractive to a fluid, while in other communities, they are reserved exclusively for water. In the latter case, the terms 'solvophobic' or 'lyophobic' and 'solvophilic' or 'lyophilic' are employed for any fluid other than water. The distinction between hydrophobic/solvophobic/lyophobic and hydrophilic/solvophilic/lyophilic is usually defined by the contact angle of a drop of water/liquid in coexistence with its gas phase at a planar wall being larger/smaller than 90 degrees. If the contact angle is exactly 90 degrees, the wall is called neutral. However, this distinction is not always practiced in the literature that strictly. Not to mention that in the chemistry community, hydrophobic and hydrophilic often means reactive and non-reactive with water. For the sake of simplicity and given the generality of some themes discussed in this review, here we use the terms hydrophobic and hydrophilic with reference to any liquid, not only water, and we will distinguish a hydrophobic from a hydrophilic surface on the basis of the contact angle formed by the liquid (see Figure 1).
Theoretical and experimental research made significant progresses in understanding intrusion/extrusion, predict the pressures at which the two processes take place, how their value depends on the chemical and topological characteristics of the material, and the thermodynamic conditions, e.g. the temperature. Intrusion/extrusion is very challenging from the theoretical point of view as the process is characterized by the presence of metastable states and free energy barriers separating them. This requires special techniques, regardless a microor macroscopic description of phases is used. In this review, we summarize the recent progresses in the field, starting from a classical macroscopic framework, the confined classical nucleation theory. This fundamental framework is extended to describe more realistic systems, including nanoscopic effects treated using both continuum and atomistic techniques.
Experimental research focused on intrusion/extrusion of relatively complex materials: random mesoporous materials, often grafted with hydrophobic coatings, and crystalline nanoporous systems, such as zeolites and metalorganic frameworks, in which cavities are produced by the regular disposition of atoms in the lattice. Here we present techniques used in the experimental research, as well as technological applications of intrusion/extrusion.
Despite the intense effort, many questions remain unresolved. They are discussed at the end of this work, together with challenges in linking theoretical and experimental research to improve the synergy between these two complementary approaches and deepen our understanding of intrusion-extrusion processes in nature and technology.

Confined classical nucleation theory of extrusion and intrusion
In this section, we introduce a formal thermodynamic framework for intrusion/extrusion based on the capillary theory of multiphase systems [24,25]. Within this framework, extrusion (drying, dewetting) of liquids from porous media or textured surfaces can be seen as the formation and growth of gas/vapor bubbles within the confining cavities of the imbibed solid. Intrusion (wetting) is just the opposite process, the shrinking of a gas/ vapor bubble initially occupying the cavities of the solid. In the following, we will refer to gas and vapor bubbles without distinguishing between the two cases; the differences between gas and vapor bubbles will be discussed in Sec. 2.2 and 3.4.
It is convenient to write the capillary model of our three phase system within the grand canonical ensemble, i.e. the ensemble at constant chemical potential, volume, and temperature. The thermodynamic potential of the grand canonical ensemble is the grand potential Ω, which for a bulk system reads Ω bulk ¼ À PV, where P and V are the pressure and the volume occupied by the (bulk) system. For a system composed by several phases, the simplest form of the grand potential is the sum of bulk and interface terms: Where P x and V x are the pressure and volume, respectively, of the phase x, liquid (L) and bubble (B). Here, it is assumed that the solid does not change along intrusion/extrusion; hence, its bulk term is constant along the process and is embedded in the arbitrary constant of the thermodynamic potential. The interface term is proportional to the area A xy between two phases x and y via the surface tension γ xy > 0. The physical meaning of the quantities appearing in Eq. 1 is shown in Figure 2. The second equality arises from the Figure 2. A) Liquid intruding a pore. Bulk and interface terms of Eq. 1 are shown graphically, together with the geometrical meaning of the contact Young angle θ Y . b) Graphical representation of the conditional (constant V B ) minimization of the grand potential vs � LB , thus obtaining the optimal liquid/bubble interface � � LB . c) Sketch of the free energy/grand potential profile along intrusion. The system is initially extruded and V B is large; to intrude, the pore the liquid must cross a grand potential barrier much higher than the thermal energy, k B T, which results in a large number of unsuccessful attempts (represented by the red arrow).
conditions V L þ V B ¼ V, with V volume of the sample available to the fluids, A SB þ A SL ¼ A S , with A S the exposed area of the solid, and cos θ Y ¼ γ SB À γ SL À � =γ LB , with θ Y the Young contact angle of the liquid. An interesting consequence of Eq. [1] is that while liquid/vapor bulk coexistence occurs when P B ¼ P L , within porous systems, coexistence depends on the contact angle and solid surface exposed to the fluids: , with A pore area of the mouth of the pore, determining the A LB when the liquid is extruded. Thus, one sees that the confined equilibrium conditions depend on the chemistry of the solid and liquid, i.e. the hydrophobicity of the material cos θ Y (see Figure 1), and on the geometrical characteristics of the pore, V, A S , and A pore .
The nucleation path in a bulk liquid within the capillary theory, the classical nucleation theory [25,26] (CNT), can be written in a variational form, searching for the liquid/bubble surface � LB minimizing the grand potential Ω for a bubble of volume V B [27]. This approach can be extended to the case of extrusion, i.e. formation of a bubble within a liquid imbibing a porous solid [27,28]. In fact, � LB directly determines A LB , and the intersection between � LB and the solid surface exposed to the fluid determines A SB . Thus, for a given Figure 2) and, consistent with the CNT, the extrusion path is the sequence of surfaces � � LB associated to bubbles of (growing) volume V B . Within this framework, in which one assumes that the bubble relaxes to its interface shape � � LB at each value of its volume V B ; intrusion is the same path followed along the opposite direction.
The framework introduced above not only allows one to predict the confined coexistence conditions but also the extrusion kinetics through the relation [25,29]: Where ΔΩ y is the (grand potential) barrier to be overcome along extrusion (intrusion) and β ¼ 1=k B T is the inverse thermal energy at temperature T, with k B the Boltzmann constant ( Figure 2). In the presence of a barrier, the system must perform multiple attempts to pass in the other basin even if this has a lower free energy/grand potential. The extrusion time τ ¼ 1=k can be obtained by inverting Eq. 1: obtained by identifying maxima of the grand potential along the extrusion (intrusion) path determined with the variational approach described above. Usually, the pre-factor k 0 is estimated by the Blander-Katz [30] or Kramers equations [25]. Once more, we remark that this is the confined counterpart of the CNT, and for this reason we call it confined CNT (cCNT).
Euler equations associated with the variational problem underneath cCNT impose that � LB has a constant mean curvature and must form an angle # Y with the porous solid [28]. The determination of the cCNT extrusion path is non-trivial, and analytical solutions were found only for a few relatively simple but relevant cases: 2D pore [28] and 3D conical crevices [31]. In all other cases, one must find numerical solutions to the variational problem min This numerical approach has been used to identify extrusion path from regular textured surfaces of complex morphologies [32][33][34][35] and cylindrical pores [9,36]. Let us illustrate an application of cCNT. Figure 3 shows the extrusion path of a liquid from a 2D pore. The dashed curves represent the LB interface of minimum grand potential at different stages along the process. At variance with bulk bubble nucleation, extrusion from a confining environment is non-trivial: i) the bubble is initially formed at one of the bottom corners, then, ii) when the bubble is large enough, it changes shape, becoming a conventional meniscus bridging the two lateral walls. Hence, iii) the triple line -a point in 2D geometries -reaches and gets stuck at the top of the cavity, and iv) the bubble eventually grows outside the pore, v) depinning from the cavity mouth only when the bubble is large enough that the Gibbs criterion -that the contact angle with the horizontal surface is equal to the Young's contact angle -is met. Figure 3 shows the barrier as a function of pressure. Figure 3 allows to reconcile the mechanical and statistical mechanics perspectives of extrusion. Since the dependence of the extrusion time on the barrier is exponential, as far as the pressure is too high and the barrier is sizably above k B T, the τ is largely above the typical experimental time, from seconds to hours. τ quickly (exponentially) drops when the barrier approaches the value k B T ( Figure 3). This change happens on such a small pressure scale to be consistent with the concept of extrusion pressure, the well-defined value of the pressure at which, according to the mechanical perspective, there is no force preventing the bubble to grow and the liquid to be expelled from the cavity.
Remarkably, cCNT predicts a morphological transition along extrusion between bubbles of different morphologies. This is seen in other systems as well [9,[31][32][33][34][37][38][39][40][41][42][43][44] and is consistent with atomistic results (see below). However, the sharp transition conflicts with the dynamic point of view, requiring a smooth evolution of the liquid meniscus � LB ; this smooth path is recovered with an extension of the cCNT [45], implementing the string method for continuum systems [46]. The continuum string reveals a path analogous to the one obtained from cCNT, differing mainly in correspondence of the morphological transition ( Figure 3).
Why developing a theoretical framework for intrusion/extrusion is important? As mentioned above, cCNT allows to identify the so-called kinetic bottleneck of intrusion and extrusion, the configuration � LB corresponding to the energetic barriers ΔΩ y . This allows to interpret intrusion/extrusion experiments, see Sec. 5, and similar to chemical reactions, design structures with controlled intrusion/extrusion barriers to speed up or slow down wetting and drying at relevant thermodynamic conditions. In other words, one can use the topology of cavities as a catalyst, if desired a negative catalyst, to either and extrusion (red) barrier as a function of the pressure. Here, for the sake of generality, both the pressure and barrier are dimensionless; Δp ¼ P L À P B ð Þ=Δp max , with Δp max ¼ À 2γ LB cos # Y =l, where l is the characteristic length of the texture (see panel a); where ΔΩ y is the transition barrier, from either Wenzel to Cassie-Baxter (W → CB) or vice versa (CB → W), and V is the volume of the texture. In particular, in this panel, we show the intrusion/extrusion barrier as a function of the pressure for a system with # Y ¼ 110 � , a typical value of the so-called chemical hydrophobicity for materials of technological interest in the field of superhydrophobicity. c) extrusion time τ ext as a function of the liquid pressure P L for a cubic pore of 10 nm characteristic length and # Y ¼ 126 � . The horizontal line corresponds to τ ext ¼ 1 s. d) Configurations along the liquid extrusion from the cavity of panel a) as obtained from the continuum string method. The blue region represents the liquid domain, which changes along the process. Adapted from Refs. [21] and [38]. accelerate or to hinder intrusion/extrusion. Finally, having a theoretical framework for intrusion/extrusion helps identifying nanoscopic effects, here identified as those effects exceeding the predictions of the capillary theory. Nanoscopic effects are discussed in the following section.

Effect of dissolved gas
So far, we did not distinguish between gas and vapor in the theoretical treatment of intrusion and extrusion. While to a large extent the simpler theory introduced above works also in presence of gases, when their solubility in the liquid is very low, some effects emerge that require an ad hoc treatment. Low soluble gases, e.g. air, are often naturally dissolved in water, and they have a high probability to be expelled from the liquid bulk and to accumulate at solid surfaces [11]. Moreover, gases are sometimes exploited to enhance bubble nucleation, e.g. in the fields of energy storage nanodevices and superhydrophobic surfaces, discussed in Sec. 5. Gas-enhanced bubble nucleation is sometimes detrimental, e.g. bubble nucleation hinders nanopore sensing [47,48] and chromatography [4,49].
As far as the gas is only dissolved in the liquid and gas atoms/molecules do not reach the cavities to be intruded/extruded, Eq. 1 still holds. On the contrary, if there is gas in the cavity to be intruded, or if the bubble bringing to extrusion is made of gas, Eq. 1 might require some amendment. The general treatment in which gas diffusion through the liquid meniscus within the cavity might be the limiting step of the intrusion/extrusion process goes beyond the framework within which Eq. 1 is derived, quasi-equilibrium. Hence, here, following ref. [32], we consider two limiting cases consistent with this framework: infinitely slow gas diffusion, in which the amount of gas present in the bubble does not change along the process, and infinitely fast diffusion, where the amount of gas is determined by equilibrium relations.
In the case of fast gas diffusion, the system can be assumed to be in chemical equilibrium. In this case, Eq. 1 must be amended by adding to P B , the vapor tension of the liquid, P g , the partial pressure of the gas. This results in a new volume contribution to the grand potential of À P g V B that can lead to a consistent decrease of the nucleation pressure, even in diluted conditions where the thermodynamic coefficients of the liquid phase in Eq. 1, P L and γ SL , will be changed only slightly by the gas atoms as compared to the pure liquid phase. Calculations of density functional theory (see Sec. 3.3) have shown that even small concentrations of xenon dissolved in liquid water are sufficient to alter the equilibrium between the filled and the empty states of a hydrophobic nanopore: the overall effect of xenon is to shift the equilibrium toward the empty state of the nanopore [10].
Let us now move to the case of infinitely slow gas diffusion. Within this hypothesis, the case in which the starting condition is the liquid completely filling the cavities/textures can still be handled with Eq. 1. In fact, for an infinitely slow gas diffusion, the bubble forming along extrusion is still a pure vapor bubble. Consider now the case that the cavities/textures initially contain a vapor and gas bubble of volume V B;0 and a gas partial pressure P g;0 . Within the hypothesis of infinitely slow gas diffusion, i.e. assuming that the liquid-vapor interface is impermeable to the gas, one must take into account the isothermal expansion of the vapor-gas bubble along extrusion. Assuming the ideal gas behavior, one must add to the grand potential the volume term P g ¼ À P g;0 V B;0 ln V B V B:0 . This term helps expanding the initial vapor-gas pocket by lowering the free energy barriers for nucleation, as it was shown for a wide conical crevice [31]. In Sec. 3.4, we shortly discuss the effect of insoluble gas dissolved in water on the extrusion barrier and kinetics as obtained from atomistic simulations; here, we anticipate that these results are qualitatively consistent with the theoretical predictions discussed in this section. Equation 1 suitably describes relatively large bubbles. For very small bubbles, the continuum description underneath Eq. 1 completely breaks down, through morphometric forms [12,50] have been developed that extend the validity of this class of descriptions [51]. When V B is on the nanometer scale, however, one can still rely on continuum models of the three-phase system, but Eq. 1 is usually modified to take into account additional effects.

Nanoscopic corrections to the cCNT
Tolman length. The first phenomenon to consider concerns the change of γ LB along nucleation and growth of the bubble. In fact, it is intuitive that the strain of the liquid at the interface, the variation of the intermolecular distance, changes with the size of the bubble. Tolman [14,16,17,52] derived the leading order correction to the surface tension depending on the curvature of the bubble (1=R in the case of spherical bubbles): Here, δ is the Tolman length, which determines the dependence of the surface tension on the curvature of the bubble. For large bubbles, R≫δ, the surface tension is equivalent to that of the flat surface (R ! 1). The Tolman length of water has been a source of debate throughout the literature, with different measurements or predictions differing in both magnitude and sign. This difficulty arises partly because the Tolman length enters thermodynamic expressions in the same manner as the line tension, discussed further below. Therefore, it is often difficult to disentangle effects of curvature from those of the three-phase contact line. Despite this, the Tolman length can be estimated from computer simulations, consistently producing δ ≈ −0.1 nm for various models of water. From Eq. 1, this means that curvature increases the surface tension of nanoscale bubbles.
Recently, Menzl et al. have shown that if one takes into account the size dependence of the surface tension, CNT quantitatively models numerical results obtained from simulations for bulk nucleation (Figure 4) [53]. Here, the quantitative matching allowed to estimate the value of the Tolman length is δÀ 0:2nm. 2 Line tension. As shown in Figure 1 and 3a, a bubble or droplet in contact with a solid surface result in a three-phase contact line where the SL, SB, and LB interfaces meet. Associated with this contact line is an excess free energy per unit length, called the line tension, τ. The line tension modifies surface thermodynamics, represented as an additional contribution τλ to the grand potential of Eq. 1, where λ is the length of the contact line. Similarly, inclusion of line tension modifies Young's equation according to γ SB À γ SL ¼ γ LB cos θ Y À τ=R τ , where R τ is the radius of the circular contact line of a spherical droplet.
The line tension has also been notoriously difficult to determine experimentally, with estimates varying by orders of magnitude, 10 −11 − 10 −6 N, as well as varying in sign [18,[54][55][56]. Others have questioned its importance, finding that quantitative agreement between atomistic simulation results and capillary theory can be obtained without including line tension effects [40,41]. Measurements of drying pressure for water evaporation in hydrophobic confinement yielded consistent values on the order of −30 pN [57]. A mechanical route for computing line tensions from molecular simulations was recently developed and predicted a line tension for water in hydrophobic confinement of −2.4 pN, a smaller in magnitude but still negative value [58]. Thermodynamic fitting of simulation results has also yielded negative line tensions closer to the experimental result [44,59]. The negative sign of the line tension impacts drying by favoring the formation of bubbles, a crucial step in extrusion processes.
In mixtures, the line tension can be tuned by varying the concentration of one component [58,60]. Line tension effects can oppose those from contact angle and surface tension changes, leading to unexpected concentration dependencies on drying thermodynamics and kinetics. For example, decreasing the salt concentration of a water-in-salt electrolyte lowers its line tension, which consequently lowers the free energy barrier to extrusion [60]. This outweighs opposing effects from changes in both contact angle and surface tension [60]. In these mixtures, solutes can adsorb to the contact line. These line-active molecules or lineactants, in analogy with surfactants, can further decrease the line tension. For example, it was recently shown that carbon dioxide preferentially adsorbs to the contact line and can decrease the line tension to −7.3 pN for a hydrophobic surface and −16.9 pN at a hydrophilic surface and consequently alter the rate of bubble nucleation [58]. Such investigations into the role of lineactants in controlling surface thermodynamics and nucleation kinetics are critical to leveraging solute effects, including gases (Secs. 2.2 and 3.4), for extrusion and intrusion in device applications.
Concerning the sign of the line tension, a final remark is in order; one imagines that a negative value of τ is unphysical, as in this case Ω would be unbound from below, and bubbles with an infinite triple line of, e.g. fractal shape, will be favored. However, if one considers that the minimum physical length scale of the continuum model must be longer than the atomistic molecular size, several Ä, an infinite fractal triple line goes beyond the physical meaning of the model. In addition, though often neglected in the simpler models, any bending of the triple line is also accompanied by bending of the liquid-vapor surface nearby, which increases the liquid/ bubble area, hence the overall free energy.

Theoretical background
Using atomistic simulations, one can compute the free energy profile of intrusion/extrusion and, with suitable extensions, include effects beyond the quasi-static picture of the thermodynamic framework introduced in Sec. 2.
In statistical mechanics, any thermodynamic potential is nothing else than the probability to observe a (macroscopic) state reported on a logarithmic scale and in k B T units (k B , in the case of entropy): is the probability to find a system in any microscopic state Γ having a value � � of the observable � Γ ð Þ. Here, δ � ð Þ is the Dirac delta function selecting points in the phase space Γ consistent with the condition � Γ ð Þ ¼ � � . m Γ ð Þ is the microscopic distribution defining the nature (e.g. Gibbs vs Helmholtz free energy) of the thermodynamic potential A � � � ð Þ associated to P � � � ð Þ. Here it is worth to remark an important difference between the atomistic and macroscopic model of the thermodynamic potential, which poses some theoretical challenges in their comparison. The continuum sets aprioristically the relation with the bubbles geometrical characteristics; for example, in the case of bulk systems, the grand potential has the following depen- is the area of a sphere in terms of its volume. If one expresses the grand potential in terms of the radius R of the bubble, one simply replaces V B with 4=3πR 3 . This is different with A � � � ð Þ as one has to apply the rule of the probability of a transformed variable, which introduces the Jacobian between two sets of variables, the derivatives of V B over R is the present example. While this transformation is simple in the case of bubble nucleation in the liquid bulk, it is non-trivial for intrusion/extrusion in complex confining media.
With reference to cCNT introduced above, the relevant observable � Γ ð Þ is the bubble volume V B . V B , however, has no trivial definition in terms of the phase space Γ. The first step to achieve this goal is to notice that for a twophase, liquid/vapor system with a sharp interface, the average fluid number density within a cavity, number of fluid particles per unit volume, is where ρ L and ρ B are the liquid and bubble bulk number densities, respectively.
ð Þ is the current number of the fluid molecules in the cavity, which depends on their atomic positions r, and N L and N B would be the number of molecules if the cavity is completely filled of liquid or vapor, respectively. Thus, within the sharp interface approximation, measuring the liquid and bubble volumes at an atomistic level amounts to count overall number N r ð Þ of fluid molecules in the cavity in the present configuration r and compare this against liquid-and vapor-filled cavity results.

Simulation techniques
In principle, one could run long molecular dynamics (MD) or Monte Carlo (MC) simulations and count the number of configurations with a number N r ð Þ of fluid molecule in the cavity and build a histogram approximating the probability density of Eq. 3. However, as shortly discussed in Sec. 2, often in relevant thermodynamic conditions, the intruded and extruded states are separated by sizable energetic barriers, which prevents the system to visit both states on the timescale of simulations, maximum several hundreds of nanoseconds. This is the well-known problem of rare events [61][62][63][64][65][66]: the transition from one state of the system to another happens once in a blue moon [61,62], and one must develop some special techniques to accelerate atomistic exploration of the state space that still allows to reconstruct the statistical properties of the unbiased system. Several techniques have been introduced to address this problem, and here we discuss two tightly connected techniques, restrained molecular dynamics (RMD) and umbrella sampling (US), which have been often used to investigate intrusion/extrusion.

Retrained molecular dynamics
The key observation of RMD is that in the definition of P � � � ð Þ, the Dirac delta function can be replaced with a smooth approximant, e.g.
tion for the given ensemble, e.g. the Boltzmann weight for the canonical ensemble. Within this approximation, the so-called mean force À Þ computed along a MD driven by the augmented potential V r ð Þ þ k B T=2σ 2 N r ð Þ À N � ð Þ 2 [27,63,67,68], where V r ð Þ is the physical potential and k B T=2σ 2 N r ð Þ À N � ð Þ 2 is the restraining term, which forces the system to visit configurations consistent with N r ð Þ ¼ N � : if the system drifts away from configurations consistent with this condition, it experiences a restoring force. Through this dynamics, one can compute the mean force À dA � � � ð Þ=d� � for a set of values N � spanning the domain between the fully extruded and the fully intruded system. Then, by numerical integration via the trapezoidal rule, one obtains the free energy profile.
As it will been shown below, in most of the cases, N r ð Þ is not sufficient to characterize the intrusion/extrusion process. Using a different language, N r ð Þ is not a suitable reaction coordinate [69,70]. This is because N r ð Þ does not distinguish between bubbles of same volume but different morphologies, i.e. one can have bubbles with different � LB but same N r ð Þ -or V B in the macroscopic theories. A natural extension is to use the (number) density field ρ r; x ð Þ, i.e. the number of fluid molecules within each of the cells centered around point x discretizing the space. ρ r; x ð Þ is a convenient � LB proxy in the atomistic description of multi-phase systems [38]. The intrusion/extrusion path and the associated free energy profile for field observables can be obtained by the string method in collective variable [71]. This approach has been successfully applied to intrusion/extrusion of simple porous systems revealing the analogies and differences with respect to the simpler description in terms of N r ð Þ (V B ) ( Figure 5) [27,40,41,45,72].

Umbrella sampling
As in RMD, US provides a way to compute the probability distribution P � � � ð Þ and, ultimately, the corresponding free energy, by enhancing the sampling of states that are rarely observed in a typical equilibrium simulation. To accomplish this task, US modifies the potential of the system by adding a �-dependent term that restricts sampling to desired regions of �space. This biasing potential modifies the systems Hamiltonian to Þ is the biasing potential that depends on the order parameter �. Biasing the system in this manner shifts the probability distribution of the order parameter as P b � ð Þ / P � ð Þe À βV bias � ð Þ , such that the original distribution is reweighted by the biasing potential. We take advantage of this reweighting in US by performing many simulations with different values of V � ð Þ to Figure 5. Comparison between the free energy profile of a single variable, number of water molecules in a square 2D cavity (Figures 2, 3), or the corresponding liquid volume, for atomistic (RMD) and continuum (CREaM) models, vs the density profile ρ x ð Þ (atomistic string), an atomistic proxy of the liquid/bubble interface � LB (interface string). Δω ¼ ΔΩ= VΔp max ð Þ, with Δp max ¼ À 2γ LB cos # Y =l, where l is the characteristic length of the texture (see Figure 3). To make the comparison between the different calculations feasible, the free energy is reported as a function of the distance from the initial state, normalized for the total length of the path (normalized arc length). Here 'distance' is understood as the Euclidean distance in the space of the variables: N and ρ x ð Þ for the atomistic calculations and V B and � LB for the continuum ones. This distance has not an immediate and intuitive physical sense; rather, it is a measure of the separation between states of the system along the process, which is more conveniently represented in its normalized form, with values contained between 0 (initial state) and 1 (final state). Adapted from Refs. 28 and 45. systematically sample the desired range of �. Then, the information contained in each biased simulation, or umbrella, can be combined using established approaches like the weight-histogram analysis method (WHAM) or the multi-state Bennett acceptance ratio (MBAR) to reconstruct the unbiased P � ð Þ over the entire range of �. In the context of drying transitions, we would like to model changes in water density. In US, as in RMD, we need an order parameter that can adequately monitor the drying process. One straightforward order parameter is the number of water molecules in the volume of interest, N v . While N v can be used to create a biasing potential, VðN v Þ, the forces resulting from a discrete variable like N v are impulsive and therefore unsuitable for use in molecular dynamics simulations. Instead, the indirect umbrella sampling (INDUS) method of Patel and coworkers prescribes a coarse-graining of N v to obtain a continuous order parameter Ñ v , in a similar spirit to the replacement of a delta-function distribution by a Gaussian distribution in RMD described above [73,74]. The coarse-grained Ñ v results in continuous forces that can be used to bias MD simulations. Moreover, Ñ v closely track N v , such that sampling the space of Ñ v also adequately samples N v , and the desired probability distribution, PðN v Þ, can be computed from the resulting umbrella sampling simulations. This INDUS approach has been used extensively to study hydrophobic solvation [75][76][77][78][79], as well as the solvation of proteins and protein complexes [76,[80][81][82][83][84]. Moreover, INDUS has been used to compute free energies of drying transitions of confined water ( Figure 6) [44,[85][86][87], ionic liquids [37,88], and water-in-salt electrolytes [60], all of which showed significant deviations from thermodynamic expectations due to the formation of vapor bubbles near the confining surfaces as N v decreases. However, for large polyatomic molecules and mixtures, N v is defined to be the number of heavy atoms or the number of a specific atom type in the region of interest and does not necessarily correspond to the total number of molecules. One could define multiple N v order parameters as well and compute multidimensional free energies. In practice, a single N v has still proved insightful for studying drying in complex liquids [37,60,79,88].

Kinetics of intrusion/extrusion, dynamical and non-equilibrium effects
The dynamics of intrusion and extrusion is of great interest but difficult to predict accurately from simulations. The simplest approach to estimating the rate constant is to compute the free energy profile and using transition state theory to estimate the rate, k, from the height of the free energy barrier, k / e À βΔA y � ð Þ , where ΔA y � ð Þ is the barrier height. However, transition state theory estimates are only accurate when � is the proper reaction coordinate. As discussed above, the reaction coordinate is difficult to determine for intrusion/ extrusion and instead order parameters like the number of molecules in a region of space or the water density are used. In these cases, dynamical effects must be taken into account in simulation estimates of rate constants. While various approaches exist for doing this, rate constants for intrusion/extrusion have been predicted using forward flux sampling (FFS) [89].
FFS is a path sampling technique that applies no bias to the system in order to accurately describe the system's dynamics. In FFS, one first chooses an order parameter that can distinguish between reactant (A) and product (B) states; for extrusion, A is the wet state and B is the dry state. Then, the order parameter is used to divide the state-space between A and B into many states or 'milestones' with different values of the order parameter. Many trajectories of a specified length are then initiated from A, and the corresponding transition probability to reach the first milestone can be computed from the fraction of trajectories that start from state A and reach the first milestone. More trajectories are initiated from the first milestone to compute the transition probability to the second milestone, and this procedure is repeated until the transition probabilities between all the milestones spanning A to B are computed. The product of all these probabilities with the flux of trajectories out of state A into the first milestone yields the rate constant for the transition from A to B. Piecewise reactive trajectories can be analyzed to identify the process' mechanism ( Figure 7). FFS has been used to compute the intrusion/extrusion mechanism in pores. Interestingly, it was found that when the system is allowed to depart from quasi-equilibrium, the intrusion/extrusion depends on the pressure [90]. Indeed, analyzing the simple case of a liquid intruding/extruding a 2D pore, one finds that in quasi-static (infinitely slow) conditions, the transition path is independent on the pressure, hence on the barrier (Figure 8), and intrusion and extrusion are symmetric. On the contrary, if one includes dynamical effects, i.e. if the transitions take place on a finite timescale like in Figure 7. Graphical illustration of the FFS method. Several trajectories are started from state A, at the center of the figure, and are continued only those reaching the first milestone, here illustrated as a prescribed distance from the center. From each note of the figure, representing the hit point of a trajectory with the next milestone, one starts several trajectories evolving according to a stochastic dynamics, which then can follow different paths. This operation is repeated until enough trajectories reach the final basin. The statistical analysis of these trajectories allows one to determine the process rate (see main text). Moreover, joining the piecewise trajectories allows one to investigate the mechanism of the process, to identify multiple "reactive" channels, if present. Adapted from Ref. 44. Figure 8. Intrusion/extrusion mechanism: comparison between the quasi-static and genuine dynamic paths as obtained from RMD and FFS, respectively. a,b) Intrusion/extrusion free energy profiles ΔΩ as a function of the number of fluid particles within the 2D pore, N, c) at a set of pressures spanning the range in which the intruded and extruded states are the stable ones, and the process to go from the less to the more stable state is either characterized by a sizable barrier or its barrierless. Pressures are reported in Lennard-Jones reduced units of measure, ε=σ 3 ; where ε and σ are the characteristic energy and length, respectively, the Lennard-Jones real experiments, the paths are highly dependent on thermodynamics conditions, and intrusion and extrusion follow different paths (Figure 8). For the 2D pore mentioned above, while in quasi-static conditions, the liquid front breaks the symmetry imposed by the pore, regardless of the pressure, with FFS intrusion taking place with the liquid advancing with a flat meniscus. On the contrary, extrusion always starts with the formation of a vapor bubble in a corner at the bottom of the cavity but the quote along at which the bubble flattens out depends on the pressure applied on the system. In other words, unexpectedly, inertia plays a crucial role also for liquids at the nanoscale.
FFS was used to compute the rate of water evaporation (extrusion) between hydrophobic walls, uncovering an exponential dependence of the rate on the distance between confining surfaces; the rate can vary by 10 orders of magnitude with a change in distance of 5 Å [29]. Moreover, FFS was also used to examine the dependence capillary evaporation (extrusion) rates on the flexibility of the confining material, in which material deformation impacts free energy in Eq. 1 [38,91,92]. Here, water was confined between model hydrophobic plates and the stiffness of the bonds between plate atoms was varied to alter the plate flexibility. It was shown that increasing plate flexibility can increase the capillary evaporation rate by many orders of magnitude -decreasing the bond stiffness by 10 kT/Å can increase the evaporation timescale from seconds to nanoseconds! As suggested by the authors, this huge response to the flexibility of the confining material could be leveraged by biological and material systems. For example, reversible changes in the chemistry or molecular conformation could alter the flexibility of a confining surface, inducing a rapid and switchable extrusion response to turn off conduction through a channel.
particles of the fluid. d) The mechanism is investigated considering the degree of symmetry of the liquid along intrusion, measured by the difference N 1 À N 2 of fluid particles in the left and right halves of the pore. Panel d) shows that regardless of the pressure, if the process is studied in quasi-equilibrium conditions (RMD, INDUS and other methods), i.e. giving the system enough time to reach equilibrium at the present level of filling N, starting from ~60% of filling, the symmetry is broken, with more liquids in the left (N 1 À N 2 > 0Þ or right (N 1 À N 2 < 0) half of the pore. This mechanism is followed regardless of the pressure, and intrusion and extrusion take place following the same path along opposite directions. (f-m) A similar analysis is shown for FFS simulations, which allow to take into account genuine dynamical effects (FFS). In this case, where we report the color map of the (logarithm) of the probability density to find the liquid in a configuration consistent with prescribed values of N 1 À N 2 for a given total number of liquid particles in the cavity, N, along intrusion (f-h) and extrusion (i-m) trajectories at various pressures. The color code of the map is reported at the bottom of the figure. One notices that in this case i) intrusion and extrusion follow significantly different paths and ii) within each process, either intrusion or extrusion, the path depends on the pressure. These results stress the relevance of dynamical effects in intrusion/extrusion and the need to go beyond the quasistatic picture. Adapted from Ref. [90].

Atomistic simulations in practice: software
Performing RMD, string, umbrella sampling, or forward flux sampling simulations require so additional software structure added atop standard molecular dynamics or Monte Carlo codes. RMD, string, and umbrella samplings need a restraining potential to be added to the physical one, associated to the collective variable(s) controlling intrusion/extrusion, � χ r i ð Þ. � χ r i ð Þ is typically a sigma-like function whose value depends on the distance of the particle i from the border of the domain. 3 This can be (and is) implemented in community codes either via 'scripting', for those codes, such as LAMMPS [93], which are provided with a powerful (input) scripting language controlling the execution of the simulations, or via external 'plug-in' codes, such as PLUMED [94][95][96] or COLVARS [97]. It is worth remarking that enhanced molecular dynamics and free energy calculations are very frequently based on much simpler collective variables, such as distance between atoms, angles, dihedrals, etc. The most used collective variables depend on a reduced list of atoms of which one computes some prescribed function(s), and the computational cost of their determination is negligible with respect to the cost of force calculation. This is not true in the case of the collective variables of our interest, e.g. N r ð Þ, as due to the diffusivity of fluids, one must compute the contribution of each particle to the collective variable. Thus, one must take advantage of all strategies typically used to compute interatomic forces, such as the linked cells [98,99]. In our experience, it resulted convenient to either use the scripting capability of codes or implement the collective variables of our interest in the code itself. 4 A similar analysis holds true for the forward flux sampling. Codes exist, which can be combined with molecular dynamics engines, e.g. the Software Suite for Advanced General Ensemble Simulations -SSAGES [100]. It is also worth mentioning the OpenPathSampling [101], a python library that can be combined with MD engines, presently GROMACS [102] and OpenMM [103]. Nevertheless, also in this case, a bottleneck is the high computational cost of the typical intrusion/extrusion collective variable(s); we found convenient to adopt our home-made implementation of the FFS within LAMMPS, obtained by a suitable combination of hacking of the original code, to efficiently implement the collective variables of our interest and a suitable control of the simulation via the facilities offered by the scripting input.

Density functional theory
Classical density functional theory (DFT) [104] can be considered a link between the macroscopic description of a confined fluid by thermodynamics and the fully microscopic account of a computer simulation, as DFT provides a functional of the grand potential with the form Ω ρ where F ρ ½ � is the functional of the intrinsic free energy, which describes the mixing of the ideal gas as well as the interparticle interaction, V ext x ð Þ is the external potential, which accounts for the confining geometry, and μ is the chemical potential. One can prove mathematically that such a functional exists and that it possesses the property that it is minimal and equal to the grand potential, i.e. the thermodynamic potential, if it is evaluated for the equilibrium density profile. This allows one to define variational principle and determine the inhomogeneous equilibrium density distribution of a system and the corresponding grand potential by minimizing the functional: where ρ 0 r ð Þ is the equilibrium density profile. This variational principle results in a Euler-Lagrange equation that has to be solved numerically on a grid.
While the starting point for DFT is, like in computer simulations, the Hamiltonian of the system, the functional F ρ ½ � typically does not represent the underlying Hamiltonian exactly but introduces some approximations. It is also instructive to note that the results obtained by DFT and by atomistic simulations do have different characters. While from instantaneous atomic configurations one obtains instantaneous densities (snapshots) of the system and has to compute time averages, in the case of MD simulations, or configuration averages, in the case of MC simulations, DFT outputs ensemble averaged densities per construction. This in turn implies that the influence of rare events on the density profile might be smeared out and difficult to identify, without the use of special techniques. Also, per construction, DFT provides the equilibrium structure ρ r ð Þ (Figure 8,9) and thermodynamics Ω ¼ Ω ρ r ð Þ ½ � on equal footing [104].
Symmetries of the confining geometry can efficiently be taken into account, and so one can reduce problems of a three-dimensional fluid close to a planar, spherical, or cylindrical wall into an effectively onedimensional calculation. However, it is possible to treat geometries with little or no symmetry in an effective two-dimensional or in a full three-dimensional geometry. Typically, a DFT calculation of a full threedimensional problem runs faster than a computer simulation but requires more computer memory, if the system is large.
DFT has been applied successfully in various problems of interest in the present context, including fluid adsorption at a single planar and structureless wall [105], to capillary phase transitions between two parallel walls, in simple [106] and in complex pores [107][108][109]. DFT provides a powerful tool to study the structure of liquid-gas interfaces [104,[110][111][112] and to gain insights in the bubble formation in highly confined system at different scales, from mesoscopic pores [51] to hydrophobic gates of biological ion channels [10,13].

Effect of insoluble gases
The theoretical analysis discussed above did not take into account some peculiarities arising by the presence of insoluble or low soluble gases. Moreover, the atomistic simulations discussed so far focus on the case of liquid/vapor systems. RMD was applied to go beyond this picture, to study gas-induced bubble nucleation inside hydrophobic nanopores [113]. Nanopores with diameter of 14 Å were immersed in water containing a single atom of argon in the geometry shown in Figure 10. Pores of this size show intruded/extruded metastabilities when immersed in pure water, without any dissolved gas [114]. Two free liquid-vapor interfaces at both ends of the simulation box along the pore axis (not shown) were left to study the mechanism of drying at constant (bulk coexistence) pressure [114]. The free-energy landscape of drying of the nanopore was computed as a function of two variables, one controlling the number of water molecules inside the nanopore, N, and the other determining the (projected) position of the argon atom along the pore axis. This landscape indicates the presence of four main metastable states ( Figure 10): wet pore and argon in liquid bulk (A), dry pore argon in liquid bulk (B), wet state and argon inside the nanopore (C), and dry state and argon inside the nanopore (D). The D state is the stable minimum for a hydrophobic pore with a Young's contact angle of Θ Y ¼ 104 � . This indicates that the presence of the gas atom inside the hydrophobic nanopore further stabilizes the dry state. The drying mechanism can follow two paths joining the state A to the state D. One connects the states A-B-D, and the corresponding free energy profile is shown in red in the right panel of Figure 10. The first step of this mechanism consists of drying of the nanopore in absence of argon (AB), and its kinetics is characterized by a free-energy drying barrier of the nanopore. The second step consists of the infiltration of the argon atom inside a dry nanopore (BD), which is a barrierless event. The second path, A-C-D, shown in black in the same figure, consists of the first step corresponding to the infiltration of the argon atom inside a wet nanopore (AC), which is still a barrierless event. The second step consists of the drying of the nanopore in presence of the argon atom (C`D), which is a mechanism characterized by a vanishingly small drying barrier: one single atom of argon is able therefore to destabilize the wet state of the nanopore by suppressing the drying free-energy barrier. The second path, characterized by two negligible barriers, is hence the most likely mechanism of the nanopore drying in presence of a low soluble gas. immersed in water (blue) containing one argon atom (red). Right: free-energy profiles along two different paths connecting state A to state D. As customary in free energy calculations, the free energy is reported as the normalized arc length (see Figure 5), with values contained between 0 (initial state) and 1 (final state). The path A-C-D is the most likely drying path, characterized by two negligible free energy barriers. The right panel is reprinted with permission from Ref. 12. The 'catalyzing' effect of gas consists in decreasing the drying free-energy barriers was also observed for a nanopore that is approximately neutral to wetting. The observed gas-induced drying is also in principle independent of the concentration of the gas in the bulk, provided that one atom of gas reaches the cavity of nanopores of comparable sizes and hydrophobicity, i.e. only the time scale of the diffusion of the gas toward the cavity depends on the concentration. Overall, this study shows that dissolved gas changes the kinetics of drying [10].

Experimental techniques to investigate intrusion/extrusion
Parallel to the theoretical and computational progresses, an intense experimental activity has been devoted to identifying several fundamental aspects of intrusion/extrusion, which put in evidence potential technological applications of these processes. These activities are described in detail in the following.

Liquid porosimetry
Liquid porosimetry is a simple and reliable technique where a non-wetting liquid is forced to enter a porous material by applying hydrostatic pressure. In this technique, intrusion pressure at which a non-wetting liquid wets the pore and the corresponding volume variation (reduction) of the overall system, liquid plus suspended porous solid, are recorded simultaneously Figure 11.
Upon decompression, drying of pores occur at extrusion pressure and the system recovers its initial overall volume due to non-wetting liquid expulsion from the pore (Figure 11). It should be noted that while intrusion is guaranteed as long as porous materials is stable under the applied pressure and pore opening is big enough for liquid molecule/atom to enter, extrusion is a more complex process in which non-wetting liquid can remain trapped inside the pores of the material in a metastable state even at atmospheric pressure ( Figure 11) [35,115]. Theory and simulations [35] predicted the existence of this metastable state due to the presence of a sizable extrusion barrier at the lowest pressure achievable in normal condition, ambient pressure. A well-known example of liquid porosimetry is the classical mercury intrusion-extrusion porosimetry, the method proposed by Edward Washburn more than a century ago [116]. A schematic representation of the intrusion-extrusion porosimeter is depicted in Figure 12.
Despite the simplicity of this technique, it provides valuable information by recording most reliable quantities, such as pressure and volume variation. In particular, intrusion pressure can be used to quantify energetic barrier, which non-wetting liquid must overcome to wet the pore and can be directly compared against theoretical predictions. In turn, this represents a validation of design criteria established on the basis of simulations to tune intrusion/extrusion characteristics as well as more fundamental aspect of the theories introduced in Sec. 2. The mismatch between extrusion and intrusion pressures, for example, is helpful at studying metastabilities involved in the microscopic mechanism of the corresponding processes [68]. Moreover, if the pore is large enough and macroscopic description is appropriate, intrusion pressure can be used to study the effective advancing contact angel thorough the Washburn equation [116]. Figure 12. Sketch of an intrusion-extrusion porosimeter. A sample of lyophobic porous solid powder is placed in a liquid (in blue) and the pumps control the pressure of a pressuretransmitting fluid (in yellow). A manometer records the pressure P during the experiment, and a transducer records changes in volume of the system. The total change in the volume is then corrected for the compressibility of the non-wetting liquid and pressure transmitting fluid, in order to obtain the volume of non-wetting liquid intruded into the porous solid [117].
From the technological perspective, intrusion/extrusion pressure defines the operational conditions of devices based on this process, as discussed in Sec. 5. For these reasons, liquid porosimetry was widely applied to 'water + hydrophobic porous material' systems.
A more advanced version of liquid porosimetry implies temperature control of the process [118,119], which introduces additional parameter into the analysis and allows exploration of the temperature dependence of the mentioned above aspects. Considering that both intrusion and extrusion pressure are temperature-dependent, the forth-and-back processes can be stimulated by temperature variation instead of pressure, which might have important technological uses as described below.

Calorimetry
Intrusion-extrusion process is accompanied by heat generation [119][120][121][122][123][124] related with breaking-formation of intermolecular bonds (or solid-liquid interface development-reduction if pore dimensions are large enough for macroscopic terminology to be used). This heat can be recorded when the setup for liquid porosimetry mentioned above is coupled with calorimetry and can provide useful information on the energetics of the process. Moreover, this heat can be potentially used for thermal energy storage applications [124]. In practice, the setup for liquid porosimetry described above can be coupled to a differential calorimeter [119,124]. Typically, differential calorimetry scheme is used for this purpose Figure 13. Such a scheme allows recording heat flux generated during intrusion/extrusion, along with temperature, pressure, and volume control/scanning. For this reason, this technique is referred to as PVT-calorimetry or scanning transitiometry [125], and it allows direct and simultaneous recording of both thermic and caloric equations of states.
Since intrusion/extrusion is barrier controlled, it is non-trivial to obtain simulation results against which to compare experimental data. Indeed, one needs to compute the heat flux along the biased simulations discussed in Sec. 3.2. Thus far, this problem has not been addressed.

Electrification
It was recently demonstrated that intrusion-extrusion process is accompanied by pronounced charge generation due to triboelectrification [6,126], i.e. due to a charge transfer between the porous solid and liquid. The resulting current can be recorded if intrusion-extrusion setup is equipped with electrodes and corresponding electrical circuit Figure 14. A constant voltage is applied between electrodes, so once a porous particle becomes charged during the intrusion/extrusion process, it is attracted to the oppositely charged electrode and electric energy can be harvested. [6,126] Indeed, the underneath mechanism is still unclear. In particular, it is not clear if the mechanism interpreting so-called (planar solid-liquid) triboelectric nanogenerators (TENGs), in which electric current is the result of the combination of two phenomena: solid-liquid triboelectrification and electric induction. More research is needed to clarify this mechanism in the case of triboelectrification in porous material.

In operando crystallography
Considering the scale of the underlying phenomenon, experimental microscopic insights into the intrusion-extrusion mechanism are highly valuable. Some of the porous media used for intrusion/extrusion are crystalline porous materials, in which the internal cavities originate from the atomic positions within the crystalline structure of the material. In these cases, in operando synchrotron and neutron scattering experiments, which allow to follow the structure of crystalline material as well as structured water inside nanopores, are subjected to diffraction experiments during the intrusionextrusion process. Such experiments provide an understanding of intruded species sights in microporous materials as well as structural response of these materials to the intrusion-extrusion process.
For example, synchrotron experiments being applied to the intrusionextrusion of aqueous electrolyte solutions (MgCl 2 , NaCl, NaBr, and CaCl 2 ) into-from several zeolites (ferrierite, chabazite) revealed that both ions and water were intruded in the pores with the concentration of the intruded solution considerably higher than the initial electrolyte solutions [127,128]. In the in-operando neutron scattering experiments, effect of flexibility [126,129] and anomalous framework response [130,131] were reviled for water intrusion-extrusion into-from several MOFs (ZIF-8, Cu 2 (tebpz)). Mentioned above experiments are currently very rare as they typically require not only large facilities such as synchrotron [127,128] and neutron scattering [126,[129][130][131] but also specialized high-pressure wet measuring cells.

Technological applications of the intrusion/extrusion process
The intrusion-extrusion process is typically hysteretic, which is manifested in a difference between pressure at which a non-wetting liquid intrudes (wets) the nanopore upon compression and pressure at which it extrudes (dries) the nanopore upon decompression. This difference defines technological applicability of intrusion-extrusion process. In particular, if the extrusion pressure Figure 14. Different configurations of a high-pressure cell for recording electrification effects during intrusion-extrusion: a) hermetic flexible Teflon capsule containing two electrodes filled with a porous material and a (non-wetting) liquid. Porous material is maintained between the electrodes via Teflon tape. Rigid Teflon spacers are introduced between electrodes to ensure constant distance between them. b) Flexible Teflon tube is filled with the porous material and the (non-wetting) liquid. The two electrodes also serve as plugs of both ends of the tube. In both configurations, a and b, the capsule/tube is immersed in a pressure transmitting fluid. Using high-pressure through, wires connect electrodes to the recording equipment [6,126].
is considerably lower than the intrusion pressure, the system behaves as a bumper or as a shock-absorber [132,133] and can be used for energy dissipation applications Figure 11. If the extrusion pressure is similar to the intrusion pressure, the system behaves like a spring (sometimes called Molecular Spring [120,123,134]) and can be used for energy storage [123], thermal actuation [129], negative thermal expansion [135][136][137], and negative compressibility [130,131] Figure 11. In the following sections, we look into each of those application niches in more detail.

Shock-absorbers and bumpers
When the extrusion pressure is lower than the intrusion pressure, it simply means that amount of mechanical energy required to compress the system is higher than the amount of mechanical energy released upon its decompression. Hence, impact, shock, or vibration, which triggers the intrusion, is being mitigated within each compression-decompression cycle. The energy dissipated per cycle is equal to the area of the hysteresis loop in the compression-decompression PV-isotherm of the system (Figure 11).
An extensive work was done by the group of Eroshenko [133,149] and the group of Suciu [147,148] to construct and test car shock absorbers based on water intrusion-extrusion into-from grafted silica ( Figure 15).

Molecular springs
There are three mechanical energy storage technologies available today: flywheel, compressed air, and pumped hydro energy [150]. The flywheel, which is certainly the most advanced technology, provides very high energy density, efficiency, and good durability [150,151]. However, such excellent performance comes with a high price. Compressed air has proven to be efficient for a large-scale energy storage, particularly when coupled with wind farms. However, availability of large storage volumes (as former mines) is required for this technology, significantly limiting its applicability [151,152]. Pumped hydro energy is a classical method for mechanical energy storage, which requires appropriate geographical location with pronounced elevation difference together with water reservoir [151,152]. All three of the abovementioned technologies are only feasible at large scale, for example, when coupled with power plants. Due to a clear lack of technological diversity [150][151][152], alternative methods for mechanical energy storage are highly desirable.
When intrusion-extrusion process is non-hysteretic (or hysteresis is small enough), it can be used for mechanical energy storage as a molecular spring ( Figure 11). Molecular springs pose several advantages compared to the abovementioned technologies. First, it can be easily scaled down or up. It does not have geographical constraints and can be used both in a mobile or stationary manner. Charging/discharging happens at constant pressure (Figure 11), which permits a constant power mode, essential for most of the applications (electricity generation in particular).
It should be noted that molecular spring behaviour have only been reported for handful of microporous materials such as zeolites [153][154][155] and recently MOFs [123] when combined with water/aqueous solutions. Further progress in this topic relies on new highly stable materials with pronounced porosity as well as on new insights into intrusion-extrusion hysteresis control. Atomistic modeling can be very helpful in this context as it might allow to develop design principles on which to base novel microporous materials.

Intrusion-extrusion triboelectric generators
Triboelectrification is a fascinating phenomenon of charge generation at the interface between two mutually displaced surfaces. It has been dreamed to be exploited for practical purposes for thousands of years. However, the Figure 15. HLS car shock absorbers: a) prototypes of shock absorbers together with the PV intrusion/extrusion isotherm of the material used to dissipate mechanical energy [133]. The prototype in lab (b) [133] and field testing (c, d) [147] with and without Caron hydropneumatic suspensions (CHS)). In panel e), it is shown a comparison between the structural characteristics of conventional oil-based and HLS-based shock-absorber [148].
complexity of the underlying physical mechanism and the corresponding knowledge gap still remain the main obstacle. This is because charge transfer, pronounced heat generation (temperature variation), and kinetics of the process are extremely challenging to measure/control as they happen simultaneously during triboelectrification. To put things in perspective, note that the earliest known written mentioning of triboelectrification dates back to around 360 BC by Plato in his Timaeus [156]. Yet, in their recent review, Lacks and Shinbrot state that for triboelectrification even ' . . . most basic foundations remain poorly understood.' [157] Intrusion-extrusion process can provide a new perspective to tribocharging both from methodological as well as technological points of view. Recently, water intrusion-extrusion process was demonstrated to be accompanied by pronounced electrification effects, which was related to solidliquid tribocharging as dielectric water moves toward dielectric nanoporous material [6,126]. Considering that nanoporous materials used for intrusionextrusion poses high surface area of ~ (400-2000) m 2 /g, the available interface for solid-liquid tribocharging is indeed very large to generate considerable amount of electrical energy. Moreover, when intrusion is endothermic and extrusion is exothermic, the intrusion-extrusion cycle is associated with harvesting of heat from the environment [120][121][122][123][124]. Thus, HLS acts as a kind of heat pump converting ambient heat and mechanical work of intrusion-extrusion into electricity via tribocharging. This paved the way for the development of an intrusion-extrusion nanotriboelectric enabled heat pump exploiting dissipated mechanical work (vibrations, for example) to extract thermal energy from the environment producing electric current.
Intrusion-extrusion triboelectric generators are currently at a very early stage of development and require deep understanding of charge carrier and transfer, heat generation and corresponding local temperature fluctuations, effects of surface morphology, etc. to reach higher technology readiness levels. On the other hand, liquid porosimetry can provide a new methodology for studying solid-liquid tribocharging by exploiting i) pronounced interfacial area of well-characterized nanoporous materials, ii) their topological variability, and iii) controlled heat of intrusion-extrusion.
The first concept of intrusion-extrusion thermal actuator can be traced back to the early works of Eroshenko using low-temperature metal alloys as a non-wetting liquid [159,160]. In these works, intrusion of melted alloy is conducted into a mesoporous silica. Next, while maintaining the pressure above the intrusion pressure, the system is cooled down below the solidus point of the alloy, so it crystallizes inside the pores of silica and hence maintains there even once the pressure is released. This corresponds to the charged state of the actuator. Finally, when this actuator is subjected to temperature increase above the liquidus point of the alloys, it melts inside the pores permitting extrusion (discharge) of a non-wetting alloy releasing the mechanical energy W ext ¼ P ext � V pores . This concept was adopted as a safety actuator at nuclear power plants [161], releasing a cooling agent/ opening an emergency valve as the temperature rises to a critical value, which is equal to the melting point of the used alloy.
Later on, it was demonstrated that thermal actuation can be realized using water and grafted silica [135]. In this case, first intrusion is achieved by applying pressure. Next, upon decompression, water maintains inside the pores in a metastable state even at atmospheric pressure (charged state), which is more convenient than crystallization to generate the charged state, making it possible to use liquids such as water [135] and aqueous solutions [162]. Finally, temperature increase drives the system out of the metastable state and provokes extrusion, releasing the mechanical energy (discharge) [135].
Finally, most recently, temperature-driven intrusion-extrusion of water into flexible MOF was explored for thermal actuation [129]. Pronounced flexibility of Cu 2 (tebpz) MOF resulted in an unusual temperature dependence of intrusion/extrusion pressure ( Figure 16). This was exploited to realize an isobaric cycle, where intrusion is triggered by cooling and extrusion is provoked by heating ( Figure 16). Such cycle has several benefits compared to the previous concepts of thermal actuation: i) the heat required to provoke intrusion-extrusion is considerably lower due to subtraction of interfacial heat of intrusion-extrusion from the sensible heat required to reduce-increase temperature of the system; ii) the output work is higher due to superposition of work of intrusion-extrusion and the work of isobaric expansion; iii) temperature difference required to realize the cycle is small, ∆T ~ 15°C, due to the sensitivity of the MOF's compressibility to temperature.

Negative compressibility
Negative compressibility (NC) is the property of materials to elongate (shorten) in at least one dimension l upon compression (decompression) [163,164]. Materials having NC are scarce and are valuable for such applications as sensors, actuators, composites with effective zero compressibility, artificial muscles, smart body armours [165,166], and acoustic shielding [167].
In a recent study, it was demonstrated that ZIF-8 Metal-Organic Framework (MOF) has a gigantic linear κ l ¼ À 1 l � @l @P and volumetric κ ¼ À 1 V � @V @PÀ 10 3 TPa À 1 negative compressibility, at least an order of magnitude greater compared to the state-of-the-art [130]. As can be seen, intrusion/ extrusion process ( Figure 17) is accompanied by pronounced expansion/ contraction of the unit cell of ZIF-8 (7a) upon compression/decompression. The mechanism behind such an unusual phenomenon was related to the presence of liquid in the windows connecting the ZIF-8 cavities, which induces a swinging of the imidazolate linkers, as was demonstrated by atomistic simulation (Figure 17). This behavior is not unique of ZIF8 but has been observed also upon water intrusion-extrusion into-from another MOF -Cu 2 (tebpz) [131], which suggests that a general strategy can be designed to obtain large negative linear and volumetric compressibilities. More flexible materials and general fundamental grounds should be explored to uncover this.

Negative thermal expansion
Similarly, to compressibility described in the previous section, negative thermal expansion (NTE) materials, for which the coefficient of thermal expansion α ¼ 1 V @V @T À � P is negative, i.e. shrink upon temperature increases, are rare but attractive from the practical point of view. Among the highest known NTE coefficient values are À 10 À 5 K À 1 for ScF 3 [168], À 9 � 10 À 5 K À 1 for Ca 0.8 La 0.2 Fe 2 As 2 [169], and À 1:15 � 10 À 4 K À 1 for reduced layered (d − f) schemes of framework response to compression before and during water intrusion, respectively. Here, the blue color represents water, while the black color is used to indicate the dimensions of the framework. (g) Theoretical evolution of lattice parameter a with pressure and (h) stroboscopic trajectory of two connected ZnIm4 tetrahedra, i.e. two tetrahedra sharing an Im ligand. One observes a counter-rotation of the two tetrahedra that makes the Zn−Im−Zn triad more colinear, which increases the Zn−Zn distance, hence the lattice size. Note, due to cubic symmetry of ZIF-8, all the lattice parameters are equal and evolve with pressure in a similar way [130].
Such mechanism was demonstrated using grafted silica [136], metalorganic framework [136], and zeolite [137], suggesting a general nature of the phenomena. Moreover, as can be seen in Figure 19 [135], thermodynamic analysis suggests that pronounced thermal expansion appears when P int;ext T ð Þ approaches the pressure of the system, while the sign of this expansion, negative or positive, depends in the sign of dP int;ext =dT; which allows controlling thermal expansion coefficient in a very broad range of both negative and positive values.

Chromatography
Reversed liquid phase chromatography (RPLC) is one of the most used techniques in analytical chemistry and represented for long time as one of the first and foremost applications of hydrophobized nanoporous materials. RPLC constitutes by far the most common variety of High Performance  Liquid Chromatography (HPLC), being especially prevalent for smaller non or moderately polar analytes [171]. The 'Reverse' designation specifically alludes to the usage of hydrophobic stationary (solid) phase in conjunction with an aqueous mobile phase carrying the analytes, whereas ordinary HPLC would imply a hydrophilic stationary phase. The most common constituents of stationary phases are constituted by alkyl-bonded mesoporous silica gels variously functionalized for specific needs (most commonly octyl-or octadecyl-silanes are employed).
At a nanometric scale, the stationary phase is constituted by a very large number of small irregular pores with an average diameter of 60-100 Angstroms (for smaller analytes) and variously interconnected. The core of a RPLC apparatus ( Figure 20) is represented by the column, i.e. the device physically hosting the stationary phase and through which the mobile phase is eluted at high pressures (commonly up to dozens of MPa). As customary in chromatography, differences among the analytes in terms of affinity with the stationary and mobile phases result in different retention times within the column and consequent separation of the analytes. More, in particular, given the hydrophobic character of the materials constituting the stationary phase, hydrophobic and less polar analytes will have a higher affinity for it, resulting in higher retention times, with hydrophilic and polar species exiting faster from the chromatography column.
Aqueous mobile phases can be typically employed in order to increase retention times of analytes, yet under certain circumstances, a sudden loss of retention is sometimes experienced. This phenomenon is apparent in experimental conditions in which the flow through the column is stopped for extended periods of time and then resumed. Recently, it Figure 20. Scheme of an RPLC chromatography. A liquid solution of the species to be analyzed is inserted, together with the eluent, in the RPLC column, which contains the stationary phase responsible for the differential retention of the various analytes. The typical output of an RPLC chromatography is a chromatogram showing the traces, peaks, obtained by the detector (e.g. a mass spectrometer) for various analytes. The area underneath these peaks is proportional to the concentration of the corresponding analyte in the original solution. The analyte is characterized by a specific retention time, the time from the beginning of the measurement at which the peak appears. In ideal conditions, RPLC should produce well-separated peaks in the chromatogram, resulting from a good separation of the various analytes.
was proposed that retention losses are caused by the dewetting of the narrow cavities of the alkyl-bonded silica pores [172,173], which is supported by the observation of a corresponding reduction of the measured void volume for the stationary phase. Indeed, when stopping the flow, the system experiences a sudden decrease of pressure from high, typically to ~100 MPa, to ambient pressure, which brings to the formation of vapor bubbles in the narrowest pores, quickly extending to large sections of the nanoporous silica gels, suddenly abating the retention performances of the analytical apparatus. Repressurization allows a (perhaps partial) wetting of the pores, with the ensuing recovery of performance. Nevertheless, in most cases, repressurization is insufficient at reestablishing initial wetting conditions and performance of the chromatographic column, probably due to the intruded/extruded metastabilities discussed extensively in the theoretical section. An alternative strategy consists in rewetting the pores with the addition of significant amounts of organic solvents to the mobile phase, which, however, is not free of drawbacks, including environmental concerns associated with the disposal of RPLC waste [174].
Understanding how dewetting and rewetting depends on, e.g. pore size, surface chemistry, temperature, and many parameters, is crucial to design novel materials for the stationary phase or strategy to prevent/ recovery dewetting. For example, Walter et al. [172] noticed how larger chromatographic retention losses were correlated to smaller pore diameters of the stationary phase, while Gritti and coworkers [4,49] focused on the effect of temperature and other parameters on the dewetting kinetics. As expected being vapour nucleation a free energy barriercontrolled process, temperature is found to dramatically speed up the kinetics of dewetting. This is consistent with experimentally observed increase of water extrusion pressure with temperature [118]. Moreover, stationary phases characterized by narrower pore sizes were found to induce faster dewetting with respect to wider pores, which is consistent with theoretical models [28,44]. Gritti and coworkers also found that air, N 2 , and O 2 enhance dewetting of the stationary phases, consistent with the results discussed in Sec 2.2 and 3.4. They also found that ordinary degassing protocols have an unsatisfactory efficiency, which can be explained within the theoretical framework recently proposed by Tortora et al. [11] on the effect of surfaces in attracting low soluble gases. Finally, complex geometries of the internal cavities of the porous solid phase, with nanometric and subnanometric interconnections among larger cavities, can enhance their hydrophobicity and, in turn, the tendency to form and trap vapor and gas bubbles, as recently shown by Amabili et al. [35].
In conclusion, understanding intrusion/extrusion and formation of bubbles in confined liquid is crucial to develop novel RPLC materials and establish efficient operative protocols tailored for specific applications and conditions. 5.1.7.1. Smart pressure transmitting fluid. For the majority of investigated systems, intrusion-extrusion hysteresis is not significantly affected by compression-decompression rate [133,139]. This property is important for energy dissipation devices, like shock absorbers, where operation at high frequencies is very common. For example, intrusion-extrusion of water into-from mesoporous grafted silica was conducted at frequencies of up to 22 Hz, without noticeable effect on the hysteresis loop [133].

Other applications
The situation is drastically different when highly flexible porous materials are used: hysteresis strongly depend on the compression-decompression rate for flexible porous material [126,175]. For rigid silica, for example, water intrusion-extrusion hysteresis H is barely changing with the rate R ( Figure 21). For the moderately flexible ZIF-8, the difference is clearly evident (Figure 21), while highly non-hysteretic molecular spring based on flexible Cu 2 (tebpz) MOF effectively turns to a shock absorber as R increases ( Figure 21). These results establish a correlation between the flexibility of the porous material and the dependence of hysteresis on intrusion-extrusion rate ( Figure 21). Presently, there is a very limited number of flexible superhydrophobic porous material suitable for intrusion-extrusion. Hence, further validation of the correlation suggested by Figure 21 is needed. Materials with a high hysteresis dependence on intrusion-extrusion rate, like Cu 2 (tebpz) MOF + water system, can be used as a smart pressure transmitting fluid: at low compression-decompression rate, the fluid acts as a transmitting liquid, and at higher compression-decompression rates, it dissipates excess vibrational energy, depending on their frequency, behaving as a mechanical version of low-pass filters.

Thermal energy storage.
As mentioned above, in some cases, heat accompanies intrusion-extrusion of a liquid in a lyophobic nanoporous material [120][121][122][123][124]. Potentially, this heat can be high enough for thermal energy storage (TES) applications [124]. Indeed, thermodynamic analysis of the intrusion-extrusion process [176] suggests that combining porous materials with pronounced surface area of more than 1000 m 2 ·g −1 with nonwetting liquids having surface entropy higher than 1 mJ•m −2 •K −1 can result in heat of intrusion-extrusion of more than 1 kJ·g − [1,124]. Achieving such extraordinary values is non trivial: i) surface entropy of this magnitude is rare and is typically observed for metal alloys, which can be used only with highly stable porous materials and may require high temperature; ii) porous materials with large surface are typically microporous ones, for which thermodynamic descriptions may break down. Still, even a sizably lower thermal energy capacity of 0.1 kJ·g −1 is still attractive for TES applications [177]. These lower values are within reach, considering that heat of 0.025 kJ·g −1 has been already reported for intrusion of water into ZIF-8 MOF [119].
The energy capacity of HLS for energy storage can be improved by creating a hybrid cycle comprising a phase change material (PCM) for the non-wetting liquid and subjecting the system to melting-intrusion-extrusion-solidification cycles under isochoric conditions [124]. This is illustrated in Figure 22, which demonstrates combination of PCM with intrusion-extrusion. As thermal energy continues to accumulate, the liquid expands and the system pressurizes. When its pressure reaches P int , the liquid will enter into the pore and thus additionally store the thermal energy as the heat of intrusion. As a secondary effect, the increase of temperature reduces the surface tension of the liquid, thus lowering the intrusion pressure. The reverse process (cooling) discharges the accumulated thermal energy as a sum of the heat of extrusion, heat of fusion, and sensible heat of cooling.

Intrusion/extrusion in biological systems.
So far, we mostly focused on technological applications of intrusion/extrusion, but these processes are relevant also in biological systems, e.g. ion channels. Bubble formation, corresponding to extrusion, has been proposed as a general mechanism for the inactivation of biological channels [10]. There is mounting evidence that the (sub)nanometric confinement in the inner hydrophobic section of some class of ion channels favors evaporation of water that ultimately regulates ion conduction across cellular membranes [178][179][180]. In fact, ions require their intact hydration shell to translocate through pores [181], and the formation of bubbles can prevent their migration through ion channels even in absence of any conformational change [182]. Extrusion has also been observed before the collapse of channels [178,183], the so-called hydrophobic collapse, suggesting that this mechanism can be instrumental also in those cases in which the closure of channels is mechanical. Recent experimental [184] and simulation works [185] have shown that hydrophobic collapse can be triggered by some process, e.g. a chemical binding, in some regions of the channel. The signal is then propagated to remote section of the channels, which undergo changes of pore diameter in specific sections of channels, resulting in the liquid extrusion. This kind of signal transmission at a distance is generally referred to as allostery in the biological community.
Finally, gas-assisted formation of nanobubbles inside ion channels may also be related to the action of inert gas anesthesia [10].

Challenges for the technological applications
As can be seen from the previous sections, intrusion-extrusion process is at the heart of a very broad range of applications, requiring specific characteristics for optimum performance. However, there are certain challenges Figure 22. Thermal energy storage using melting-intrusion-extrusion-solidification cycle [124]. Initially, the PCM is solid and the porous material is extruded. Heat from the environment raises the temperature of the storage device above the melting temperature of the PCM, which melts. Any further increase of the temperature raises the pressure of the system, and when this reaches P intrusion , the porous material gets intruded by the molten non-wetting PCM. This last step increases the energy that can be stored in the system. The heat stored along this 'charging' process is released to the environment along the reverse, 'discharging' path. relevant to all of them, such as stability of porous materials, controlling intrusion-extrusion pressure, and the corresponding hysteresis and the energy balance of the process.
Stability of porous materials. Intrusion process provides challenging conditions for porous material. High pressure and forced penetration of intruding liquids often have an irreversible effect on the framework, leading to recrystallization [119,186], defects formation [121,187], or decomposition [188] of the porous structure. This limits the range of suitable porous materials, which in real applications must resist to thousands of intrusion-extrusion cycles. Grafted silica demonstrates excellent intrusionextrusion stability, which permitted reaching a high technology readiness level in the field of mechanical energy dissipation [133,147,148]. The relatively large pore size of silica (mesoporous), however, results in pronounced intrusion-extrusion hysteresis, limiting their applicability to energy dissipation. Some zeolites featuring micropores present non-hysteretic intrusionextrusion [153][154][155], which allows their employment for energy storage and conversion. However, most of them are prone to form defects upon intrusion [121,187], limiting their durability. Finally, the third and most recent class of porous materials considered for intrusion-extrusion is MOFs. For the moment, only four MOFs were demonstrated to be stable under intrusion-extrusion, namely, ZIF-8 (up to 60°C) [119], ZIF-67 (at 2°C) [188], ZIF-71 (at least at room temperature) [142], and Cu 2 (tebpz) (at least up to 90°C) [123]. Presently, the stability of MOFs upon intrusion-extrusion is still under scrutiny. In particular, some MOFs are highly stable under the action of both pressure and elevated temperature, and they tend to degrade at much lower values when intrusion-extrusion is involved. For example, ZIF-8 was reported to be stable under hydrostatic pressure of up to 0.34 GPa [189], a non-hydrostatic pressure of 1.6 GPa [190], and temperature of up to 500°C [191], but degrade at 90°C and 30 MPa when intrusion-extrusion takes place [119], including recrystallization into different allotropic forms [119,186,192]. The development and optimization of MOFs for the technological applications highlighted above requires a better understanding of the degradation mechanism upon intrusion-extrusion to improve their stability. Indeed, recrystallization upon intrusion-extrusion can be viewed as a potential mechanochemical method for synthesis of MOF polymorphs under pressures/temperatures much lower compared to non-penetrating conditions.
Understanding intrusion/extrusion pressure and their hysteresis. As already discussed, intrusion-extrusion hysteresis ( Figure 11) defines the technological application domain of HLSs: storage by molecular springs or dissipation by shock absorbers/bumpers. At the same time, the absolute values of intrusion and extrusion pressures set the operational conditions for charging and discharging. With this in mind, several strategies have been developed to tune the intrusion-extrusion hysteresis such as replacing pure water with salt solutions [142,[193][194][195], tuning the liquid viscosity [115], altering the topology of lyophobic porous materials [35], and playing with their flexibility [126]. Despite the theoretical and empirical progresses, an accurate model quantitatively relating the intrusion-extrusion pressures to the characteristics of the porous material (chemistry, pore size, topology, and flexibility) and of the non-wetting liquid is still missing. Such model would greatly contribute to advance the field, making it possible to select or design materials for specific engineering applications.
Energy balance for intrusion-extrusion process. Most of the works on intrusion-extrusion focus on the mechanical energy involved in the process, neglecting the thermal and electrical contributions to the energy balance, which, on the contrary, can be not only significant but even predominant. The heat of intrusion-extrusion highly depends in both sign and magnitude on the topology of the porous solid [196], viscosity of the liquid [115], temperature [119], and interfacial properties [124]. In ZIF-8 + water, for example, heat of intrusion was reported to be more than three times higher compared to mechanical energy [119]. Electrical energy involved in the intrusion-extrusion cycle was just recently demonstrated [6,126], and its origin, proposed to be related to liquid/solid triboelectrification, [6,126] still must be confirmed. Understanding interplay between these energies in the intrusion-extrusion cycle is essential from both technological and fundamental point of view. For example, how heat generation affects triboelectrification or how electrification affects wetting (intrusion-extrusion pressure) is essential.

Conclusions and outlook
We have shown that liquid intrusion/extrusion in and out of porous media is a process of great interest with many potential applications. Here we focused on technological applications, but this process is relevant also in biological systems, e.g. ion channels. Despite its relevance, for a long time the laws and system properties affecting the basic characteristics of the process have been poorly understood. In this article, we revised the recent research in the field, summarizing the achievements of the latest 10-15 years. Progresses on the experimental and theoretical aspects nicely complement each other, and development of stable hydrophobic nanoporous materials, allowing multiple water intrusion/extrusion cycles, permit a direct comparison between experiments and simulations.
Despite the relevant progresses, many questions remain to be address, and some of them, like liquid/solid electrification and thermal energy, call for further major progresses before a thorough theoretical description of the process will be available. Other challenges concern the chemical complexity of some of the emerging materials, such as MOFs and COFs. From the point of view of simulations, this requires the development of interaction model of the level of accuracy of ab initio model. However, the cost of this type of simulations would be prohibitive; hence, this seems the ideal playground for machine learning potentials [28,197,198]. From the experimental and applicative perspective, a crucial problem to address is to develop novel nanoporous material able to withstand cyclic intrusion/extrusion: only a handful of present materials possess this characteristic.

Notes
1. Here we used the loose notation Ω � LB ; V B ð Þ, while, more formally, we should have used square brackets denoting functionals rather than functions of � LB . 2. In the original article, the author assigned to the Tolman length a positive sign. This depends only on the sign convention of the equation for γ LB R ð Þ, which is opposite from that used in the present manuscript. 3. A form of � χ r i ð Þ converging to the characteristic function can be obtained by positioning a Gaussian function of width α, g α r À r i ð Þ ð Þ, on each fluid particle and defining � χ r i ð Þ ¼ ò domain drg α r À r i ð Þ ð Þ: in the limit of infinitely narrow Gaussian, α ! 0, � χ r i ð Þ converges to χ r i ð Þ. If α is sufficiently narrow with respect the size of the domain within which we want to count the number of fluid particles, ò domain drg α r À r i ð Þ ð Þ can be rewritten as the product of error functions of the distance of the particle i from the border of the domain. The error function is not numerically efficient to handle on a computer and is, hence, replaced by some convenient approximation, e.g. the Fermi function 1= exp λd ð Þ þ 1 ð Þ, with d the distance of the particle from the border of the domain and λ the parameter controlling the smoothness of the fermi function. 4. RMD, String, INDUS have been implemented in 'modified' versions of LAMMPS, GROMACS, other codes. Reader interested in running calculations of the kind described in this review and having insufficient coding experience are welcome to contact the authors to receive support.

Disclosure statement
No potential conflict of interest was reported by the author(s).