Bridging scales in disordered porous media by mapping molecular dynamics onto intermittent Brownian motion

Owing to their complex morphology and surface, disordered nanoporous media possess a rich diffusion landscape leading to specific transport phenomena. The unique diffusion mechanisms in such solids stem from restricted pore relocation and ill-defined surface boundaries. While diffusion fundamentals in simple geometries are well-established, fluids in complex materials challenge existing frameworks. Here, we invoke the intermittent surface/pore diffusion formalism to map molecular dynamics onto random walk in disordered media. Our hierarchical strategy allows bridging microscopic/mesoscopic dynamics with parameters obtained from simple laws. The residence and relocation times – tA, tB – are shown to derive from pore size d and temperature-rescaled surface interaction ε/kBT. tA obeys a transition state theory with a barrier ~ε/kBT and a prefactor ~10−12 s corrected for pore diameter d. tB scales with d which is rationalized through a cutoff in the relocation first passage distribution. This approach provides a formalism to predict any fluid diffusion in complex media using parameters available to simple experiments.

F luid diffusion in porous media involves complex phenomena arising from the restricted diffusivity imposed by the host porous geometry and the fluid/solid interaction [1][2][3][4] . While the medium morphology and topology impact the fluid dynamics at almost any pore lengthscale d, the effect of fluid/solid forces roughly scales with the porous surface to volume ratio S/V~d −1 5, 6 . This leads to rich dynamics in nanoporous media (for which d is of the order of the intermolecular force range ζ) with intriguing aspects such as anomalous single-file diffusion, intermittent Brownian dynamics, stop-and-go diffusion with an underlying surface residence time, etc. 7 .
For simple pore morphologies (e.g., planar, cylindrical), a unifying picture has emerged with well-identified dependence on temperature T, fluid density ρ, mean free path λ, pore size d, fluid molecule size σ, etc. 5,8,9 Single-file diffusion is restricted to d~σ while the diffusion mechanism for d ≳ σ depends on ρ and λ; Knudsen diffusion for fluids with λ ≫ d and molecular diffusion for λ ≲ d. For materials with large S/V, diffusion involves intermittent dynamics with subsequent surface adsorption and in-pore relocation steps 10,11 . When relocation is negligible (i.e. at low T and/or ρ where the pore center is depleted in fluid), diffusion is governed by surface diffusion described using the Reed-Ehrlich model 12,13 . In contrast, when relocation contributes to the overall dynamics (non-negligible pore center density), the intermittent Brownian motion is a rigorous formalism to upscale the local microscopic dynamics to any upper scale 14 . Diffusion in disordered porous media is far more complex as coupled geometrical and surface interaction effects lead to novel phenomena [15][16][17] . The fluid diffusion in such heterogeneous solids involves a non-trivial diffusivity landscape as surface diffusion/in-pore relocation boundaries are ill-defined. Diffusion in such rough landscapes is even more puzzling for nanoporous media as (1) the underlying propagatorsi.e., the probability that a molecule moves by a quantity r in a time tare not necessarily Fickian 18 , and (2) nonvanishing surface interactions in the pore leads to self-diffusivity D s different from the bulk even far from the surface 19 .
Due to the continuum hypothesis breakdown at the nanoscale 1,16 , statistical mechanics is the appropriate formalism for complex diffusion in disordered media 7 . In particular, generalization of molecular intermittence to heterogeneous media using the Fokker-Planck or path integral formalisms allows linking microscopic to macroscopic dynamics 20 . However, while these approaches rely on available material parameters (e.g. porosity ϕ, S/V ratio, structure factor S(q)), fluid dynamics concepts such as surface residence, in-pore relocation, and their time constants are often used as guessed inputs (typically, relocation/ surface diffusion are assumed to be Fickian with diffusivities equal or orders of magnitude slower than the bulk 14 ). While this qualitatively captures the complex dynamics at play, there is a strong need to establish physical laws from simple parameters such as pore size d and fluid/solid interaction strength ε. In this context, hierarchical simulations 13,21,22 allow upscaling the microscopic dynamics assessed from atom-scale simulations into kinetic Monte Carlo simulations; a precalculated free energy map ΔF(r) is used in a random walk approach with corrected hoping rates k $ exp½ÀΔFðrÞ=k B T 23,24 . However, extension to disordered solids is almost intractable because of their large representative elementary volume. Moreover, despite their robustness, such extensive simulations do not provide simple laws based on d, T, ε, ϕ, etc. because they are performed for a peculiar system under some given thermodynamic and dynamical conditions.
Here, we address the problem of fluid diffusion in ultraconfining disordered nanoporous materials by reporting robust physical laws established in the framework of surface/pore diffusion intermittence. By mapping molecular dynamics (MD) simulations onto mesoscopic random walk (RW) calculations accounting for surface residence, our hierarchical approach captures the fluid diffusion in disordered nanoporous media and their underlying complex diffusivity landscapes. Moreover, by varying the matrix porosity ϕ and pore size d but also the fluid/ solid interaction strength ε, the proposed approach provides a means to quantitatively bridge the microscopic and mesoscopic dynamics in such complex environments using simple parameters. Both the typical surface residence and relocation timest A , t Bare found to derive from physical laws involving the pore size d and the fluid/solid interaction strength normalized to the thermal energy ε/k B T. In more detail, t A is shown to obey a transition state theory t A $ t 0 A expðÀΔF=k B TÞ where ΔF~ε is the free energy barrier that must be overcome to escape from the interaction field generated by the solid and 1=t 0 A is the characteristic escape attempt frequency. t 0 A is found to be of the order of~10 −12 s (a commonly accepted value) with a correction that accounts for pore diameter d/ξ (with ξ~σ, i.e. the molecule size). As for the relocation time t B , it is shown to scale with d as quantitatively predicted by introducing a time cutoff t c~d 2 /D 0 in the relocation first passage probability distribution.

Results
Our coarse grain model is developed in the spirit of the continuous time random walk (CTRW) as first proposed by Montroll and Weiss 25 and later extended by Shlesinger and Klafter 26 to the Levy walk model and other variants. The intermittent dynamics proposed in our approach involves a waiting time distribution at the pore surface coupled to a bridge statistics taking into account the first passage probability to connect one point at the interface to another through a random walk in the accessible pore network. The latter statistics couples distance and time as in the Levy walk (coupled memory) which also pertains to the Knudsen regime 27 . In the present approach, we will mainly consider the time distribution of these bridge statistics which are mapped onto atomscale dynamics simulations to establish a bridge between the microscopic and mesoscopic scales. While the mapping proposed in this paper is derived for a simple fluid confined in prototypical models of highly disordered materials, we believe it can be extended to a much broader class of fluid/solid couples. However, such generalization must be performed with caution as there are a number of limitations which can lead to departure from the simple intermittent Brownian motion at the heart of our approach. Depending on the nature of the confined fluid and host solid, different molecular interactions are at play which are either short (e.g. dispersion, repulsion) or long (e.g. electrostatic, polarization) ranged. While such molecular interactions often lead to similar confined diffusivity mechanisms, they can induce more complex behaviors that are not entirely captured by a simple stop-and-go process. Moreover, for host solids with spatially-extended pore correlations (e.g. fractal solids), additional complexity and/or additional specific effects are expected.
Different topological porous media. Fluid diffusion in disordered nanoporous media was investigated by considering a set of 13 heterogeneous carbon structures with different densities ρ s , porosities ϕ, and pore sizes d. These structures-referred to as CS x with x the density ρ s ranging from 0.5 to 1.4 g/cm 3 -were created using a quenching procedure (see Methods for full details). Typically, using a cubic box of length 100 Å containing~25,000 to~70,000 carbon atoms depending on ρ s , molecular dynamics in the NVT ensemble was used with the reactive empirical bond order potential 28 in LAMMPS 29 to allow for bond formation/breaking during a 5 ns quench from 3000 K to 300 K. As an example, Fig. 1a shows the sample CS 0.70 filled with fluid molecules at their boiling point (as described in the Methods section, the adsorbed fluid density was estimated using standard Monte Carlo simulations in the Grand Canonical μVT ensemble). Figure 1b shows the porosity ϕ as a function of ρ s where ϕ is determined using a Monte Carlo algorithm; by inserting N probe molecules at random positions inside the simulation box containing the porous structure, the porosity can be estimated as the ratio ϕ~N v /N (where N v is the number of probe molecules that do not overlap with any of the porous structure atoms). For each structure, provided a large number N v is considered, the pore size distribution f(r) can be assessed from the diameter of the largest sphere containing each of the N v points (Supplementary Fig. 4) 30,31 . As expected, both the porosity ϕ and the mean pore size d = ∫rf(r)dr decrease upon increasing ρ s with ϕ ∈ [~0.1,~0.5] and d varying from a few to~15 Å.
Only structures with connected porosity were retained to investigate multiscale diffusion as unconnected porous samples necessarily yield zero self-diffusivities in the long time limit. For each sample, as described in the Methods section, the connectivity of the porous subspace accessible to a diffusing molecule was determined using the retraction graph associated with the digitized pore network (which conserves the topology at all scales). Such digitized binary sets are used to compute the connection number [32][33][34] : where α 0 and α 1 are the number of vertexes (either isolated or connected) and links, respectively. c t , which is a simple intensive parameter related to the number of irreducible paths per vertex, is invariant under any continuous pore network deformation. For structures with no isolated vertexes, the number of vertexes is smaller than the number of links, i.e., α 0 < α 1 , so that c t > 0 with a value that increases with pore network connectivity-in this case, the average number of links around a connected vertex is given by <N c > ¼ 2ðc t þ 1Þ. On the other hand, for poorly connected pore networks, α 0 > α 1 so that c t < 0 with c t ∈ [−1, 0]-in this case, <N c > ! 0 as c t → −1 so that the topological structure reduces to a set of isolated vertexes. In the above picture, the crossover c t = 0 is generally assumed to correspond to a percolation threshold 32,33 . Figure 1b shows that c t > 0 for ρ s ≤ 1.0 g/cm 3 as expected for a connected pore network (although c t is lower than typical values for very open networks c t~0 .5−0.7 32,33 ). On the other hand, c t < 0 for ρ s > 1.0 g/cm 3 , therefore indicating a long-range network disconnection for these dense porous structures.
Intermittent Brownian motion with underlying stop-and-go diffusion. Diffusion in the disordered media with connected porosity (c t > 0) was investigated using MD for a simple Lennard-Jones (LJ) fluid at constant temperature and for varying fluid/solid interaction strengths. For such subnanoporous materials with strongly disordered pore morphologies, provided the number of adsorbed/confined molecules is low enough, the selfdiffusivity D s is close to the collective diffusivity D c as cross-terms between fluid molecules are negligible (because fluid-solid interactions largely prevail over fluid-fluid interactions) 16 . As a result, due to the formal equivalence between permeance and collective diffusivity, i.e. K = D c /ρk B T~D s /ρk B T, the self-diffusivity also provides key insights into transport mechanisms under flow conditions as induced by pressure/chemical potential gradients. The LJ fluid parameters (σ 0 = 3.81 Å, ε 0 /k B = 148.1 K) were chosen to match those for methane-a simple nearly spherical probe. An isotropic molecular model-known as the united atom model-was used to describe the methane molecule. Such a simplified model was selected as it simply corresponds to a Lennard-Jones potential that is representative of a broad class of atomic and molecular liquids. Despite this simple fluid hypothesis, we believe that our approach can be extended to more complex fluids such as dipolar molecules. In particular, even if complex molecular structures lead to richer surface thermodynamics behavior with strong adsorption in specific sites and/or relocation with large inherent activation energies, the present approach remains relevant as such complexity is embedded-at least in an effective fashion-into the mean relocation and residence times. For each disordered porous structure, different fluid/ solid strengths were considered: ε/k B T = n with n = 0.01, 0.1, 0.2, 0.3, 0.5, 0.8, and 1. Varying ε drastically affects the porosity seen by the confined fluid since it also modifies the repulsive interaction contribution-thus inducing large changes in the effective diffusivity of the confined fluid. To probe fluid dynamics at constant porosity while scanning a broad range of ε, the fluid/ surface LJ potential was modified using a smoothing procedure involving a sigmoid function to rescale the potential well-depth at nearly constant repulsive contribution (see Methods for details).
The temperature was chosen equal to T = 450K~3ε 0 /k B to ensure that the Fickian regime is reached in all cases over the typical simulation run length (20 ns). Supplementary Fig. 6 shows the mean square displacement <jrðtÞ À rð0Þj 2 > as a function of time t for methane confined in the different disordered nanoporous materials (only data for ε/k B T = 0.1 are shown for clarity). Typically, for the disordered materials considered here, the Fickian regime is reached after a few ns as each molecule diffuses over a length scale of the order of the simulation box size L~10 nm. While such convergence is reached within typical timescales probed using molecular dynamics for these disordered materials with connected porosity (c t > 0), there are materials classes where the long-time limit extends to much longer timescales. This includes solids with longrange pore correlations such as in fractal media or strong persistence length such as in one-dimensional pores. As shown in the inset in Fig. 2a, for all systems, the self-diffusivity D swhich is inferred from the Fickian regime in the long time limit D s ¼ lim t!1 <jrðtÞ À rð0Þj 2 >=6tis lower than the bulk selfdiffusivity D 0 s . As a result, the tortuosity τ MD ¼ D 0 s =D sdefined as the ratio of the bulk to the confined self-diffusivities-is larger than 1 as shown in Fig. 2a. As expected, upon increasing ε/k B T, the average fluid/surface energy 〈U fw 〉 becomes more negative (attractive) so that τ MD increases due to the increased tortuosity adsorption/residence contribution. Moreover, τ MD increases upon increasing the solid density ρ s as more severe confinement leads to smaller diffusivity (as shown in Fig. 1, the pore size d $ ρ Àx s with x~1). The underlying microscopic diffusion mechanism in such ultra-confining materials can be identified by computing the self-correlation function G s (r, t). In particular, in an isotropic medium, 4πr 2 G s (r, t)dr is the probability distribution that a molecule moves by a distance r over a time t. As shown in Fig. 2b (see black dashed line), upon averaging over all molecules and time origins, the mean square displacement <jrðtÞ À rð0Þj 2 > displays a smooth behavior from which a confined self-diffusivity can be derived. Yet, Fig. 2b reveals that 4πr 2 G s (r, t) displays a complex behavior characteristic of stopand-go processes where the molecules switch from one location to another through jumps (the data shown here correspond to the sample CS 0.70 but the same data can be found in Supplementary  Fig. 8 for different samples and fluid/surface interaction strengths). In more detail, the probability distribution exhibits marked vertical stripes indicating that molecules tend to remain within the same spatial domain over a given time. The distance between two stripes, which roughly corresponds to the fluid molecule size σ 0 , corresponds to the jump amplitude. The typical residence time at a given position is given by the decay along the t axis. Such stop-and-go diffusion was already reported by Sahimi and coworkers 35 in molecular dynamics of gas diffusion in a carbon nanotube/polymer composite and, more recently, by Kulasinski et al. for water diffusion in amorphous hydrophilic systems 36 .
To shed more light into the complex diffusivity landscape in such disordered porous media, a single trajectory RðtÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðrðtÞ À rð0ÞÞ 2 q is provided as an example in Fig. 2b together with a visualization of the corresponding molecular trajectory in Fig. 2c (to interpret these different space-time domains, Supplementary Fig. 9 provides additional individual trajectories). Such individual trajectories are typical but not necessarily fully representative as they were chosen to identify well-defined steps. However, the mechanisms discussed below are common to all molecules and lead to the heterogeneous behavior observed in G s (r, t). Before going into details, we define here a cavity as a portion of the pore network of size d. The first narrow stripe corresponds to molecules located in a given site with displacements over short distances r much smaller than the molecule size σ 0 . Such motions are illustrated by the green portion of the individual trajectory in Fig. 2b (in this specific example, the molecule is adsorbed in the vicinity of the host surface as shown in c). This analysis is confirmed by the fact that the typical residence time associated with this dynamical sequence increases The multicolor solid line shows a typical yet representative displacement R (t) for a given molecule. Each color region (green, blue, orange, and red) indicates a portion of the 150 ps trajectory shown in (c) (see detailed description in the main text). For a more direct visualization, the corresponding trajectory is shown in (c) with successive molecular positions shown as spheres of the corresponding color. The carbon-carbon bonds in the porous structure are shown as gray sticks.
with increasing the fluid/surface interaction strength ε (Supplementary Fig. 8). The second narrow stripe centered at about r4 Å(~σ 0 ) corresponds to molecules jumping to a neighboring site. As illustrated in the individual trajectory (blue portion), such displacements correspond to molecules relocating from an adsorbed site to another while remaining within the same cavity (r < d). The third narrow stripe centered at about r ≲ d corresponds to confined diffusion where molecules explore both the pore center and surface region but remain within the same cavity (as illustrated with the orange portion of the individual trajectory). Finally, upon further increasing the time, the displacement becomes larger than the pore sizer > das the molecule is transferred from a cavity to another (as illustrated in the corresponding red portion of the individual trajectory shown in c). As expected, at large times/distances, typically when r > d, the probability distribution becomes more homogeneous as the detailed structural footprint of the disordered host matrix averages out into a single effective parameter corresponding to the tortuosity. In particular, in the long time limit, the dynamics reach a Fickian regime as the molecules diffuse over distances r large enough compared to the pore size d.
Despite the intrinsic complexity of diffusion in such rough energy landscapes, the local, i.e. in pore, self-diffusivity can be derived formally using effective approaches 19,37,38 . By in-pore diffusivity, we refer here to the short time range where molecules remain within the same cavities while reaching a pseudo-Fickian diffusion regime (in other words, such transport coefficients at the pore scale do not include network effects such as tortuosity but contain the fingerprint of the pore geometry/morphology). In more detail, considering the mean-square displacements shown in Supplementary Fig. 6, D p s can be assessed from the linear regime observed in the short time scale where <jrðt d Þ À rð0Þj 2 > ≤ d 2 (where t d is the time required to displace molecules over a distance equal to the pore size d). To further validate the inferred value, it was checked that it is consistent with the in-pore diffusivity estimated as~d 2 /6t d . The simplest effective framework consists of writing the effective pore-scale diffusivity D p s as an average over the whole pore volume, D p s ¼ 1=N R ρðrÞD s ðrÞdr where ρ(r) and D s (r) are the local density and self-diffusivity at a position r. Within the transition state theory, the bulk selfdiffusion coefficient can be written as an activated process D 0 s $ exp½ÀΔF 0 =k B T where ΔF 0 is the activation free energy to set the molecules in motion. Here, we refer to the bulk phase taken at the same temperature but also the same density as the confined phase. Therefore, even if the bulk phase is a low-density gas (for which diffusion does not involve any activation energy), D 0 s should be understood as the liquid-like diffusivity of the bulk fluid taken at the same liquid-like density. For a confined fluid, the activation energy for diffusion can be assumed to correspond to the bulk activation energy augmented by the fluid/surface potential ζ(r), ΔF = ΔF 0 − ζ(r) (the sign minus is due to the fact that the interaction potential is attractive and, hence, negative so that molecules are trapped in deeper energy sites with an escape time requiring a larger activation energy). With this assumption, D s ðrÞ ¼ D 0 s exp½ζðrÞ=k B T 19 . For complex media, there is no simple expression for ζ(r) but we use here a simple form where ζ(r) is constant when the distance to the surface is smaller than σ and decays exponentially beyond. Such a generic form leads to the following local self-diffusivity: s is the surface self-diffusivity in the vicinity of the pore surface while D 0 s is the bulk, i.e. unconfined, self-diffusivity). This expression simply assumes that the self-diffusivity is equal to the surface diffusivity D s s for distances within a critical size σ from the surface while it decays exponentially with a characteristic lengthscale r 0 towards the bulk diffusivity D 0 s as the distance to the surface increases, as depicted on Supplementary Fig. 11. After a little algebra, assuming the pore density is homogeneous, i.e., ρ(r)~ρ, one arrives at: As shown in the inset of Fig. 2a, the above effective expression provides an accurate description of the simulated in-pore diffusivity D p s . Both the variations in pore size d and fluid/surface interaction strength ε are accurately captured. The parameters D s s , D 0 s , r 0 and σ extracted from the fit against Eq. (2) can be found in Supplementary Fig. 7. As expected, the bulk self-diffusivity D 0 s is found to be constant at a value 14 ± 0.6 × 10 −9 m 2 /s. While the confined fluid density is an ill-defined quantity that depends on a given pore volume definition, we note that the bulk reduced density ρ * = ρσ 3 needed to match the bulk self-diffusivity D 0 s ¼ 14 ± 0:6 10 À9 m 2 /s inferred from this simple in-pore diffusivity model falls within the range [0.8-1] (see Supplementary Fig. 12 showing the self-diffusivity of bulk methane as a function of density at the temperature considered here). Recalling that the number of confined fluid molecules was obtained by filling each porous material at the fluid boiling point, such reduced densities further support the use of a simple effective model for the in-pore diffusivity as they correspond to typical liquid densities. Similarly, σ is independent of d and ε with a constant value of 2.4 ± 0.1 Å so that the critical distance σ for surface diffusion roughly corresponds to the fluid molecular size. Interestingly, the quality of this effective in-pore diffusivity model shows that the surface diffusion D s s and scaling r 0 can be treated as constant parameters for a given ε. On the other hand, as expected, D s s is found to decrease upon increasing ε while r 0 increases upon increasing ε. Typically, upon varying ε/k B T from 0.1 to 1.0, D s s decreases from 3 to 0.5 × 10 −9 m 2 /s while r 0 increases from 0.8 to 2 Å. The fact that the scaling parameter r 0 depends on ε can be rationalized as follows. Even if the surface/fluid interaction potential decay is independent of ε, it generates a free energy landscape ζ(r) that includes many body-fluid/fluid and fluid/wall-effects which lead to an effective scaling r 0 that depends on ε.
While the combination rule above provides a quantitative description of the molecular dynamics data, it remains mostly effective as it relies arbitrary choices combined with an empirical description of the diffusivity landscape explored by the fluid molecules. First, ζ(r) should be seen as an effective free energy field that modulates the bulk self-diffusivity by accounting for local intermolecular interactions but also for local density/ packing effects. Therefore, even with simple pore geometries, instead of a robust free energy field rigorously derived from intermolecular interactions, ζ(r) is an effective function which is used to describe the self-diffusivity decay upon increasing the distance to the pore surface. The constant surface diffusivity at the pore surface is used to account for the fact that adsorbed molecules explore homogeneously the surface region~2σ. Moreover, even if the conclusions above are qualitatively independent of the different assumptions involved, the decomposition into surface and bulk-like diffusions is also sensitive to the exact scaling defined in Eq. (2) and the parameter 2σ used to define the surface layer. In particular, other efficient decomposition rules have been proposed such as a simple weighted sum of surface and volume diffusivities which was found to accurately describe the dynamics of water in nanoconfinement 39 . Moreover, such surface/volume partition and the resulting predictions in terms of in-pore diffusivities D p s are also dependent on the geometry choice-usually far from any realistic description-made to describe the pores in such disordered materials (planar, cylindrical or spherical). In practice, as will be shown in the rest of this paper, to avoid relying on such effective frameworks, the intermittent Brownian formalism mapped onto molecular dynamics data provides a means to describe stop-and-go processes in such disordered and ultraconfining materials without invoking any definition for the surface layer and the self-diffusion decay as molecules get closer to the pore surface.
The stop-and-go, i.e. intermittent, diffusion observed in our atom-scale dynamics simulations suggests that the corresponding data can be analyzed using the framework of intermittent Brownian dynamics. Indeed, as shown in Fig. 2b, while ensemble averaging over each molecule leads to a Fickian regime with an effective self-diffusivity, each individual trajectory involves intermittent motion with alternate series of in-pore diffusion and surface adsorption. In more detail, within this formalism, the mesoscopic, i.e., coarse-grained, dynamics beyond molecular time and length scales is governed by two parameters: the residence time t A during which a molecule remains adsorbed to the surface and the relocation time t B between two adsorption periods 10,14 . To probe such intermittent dynamics, the pore space Ω available for the dynamics of spherical molecules inside the carbon matrix was extracted by mapping a 3D lattice network having a voxel size Δ = 0.2 Å (as explained in the Methods section). A voxel belongs to the pore space if its distance x to any carbon center is x > σ where σ is the LJ parameter for the fluid/ surface interaction. The surface boundary ∂Ω of Ω is made of surface voxels which are at the frontier between Ω and its complementary space. This allows defining a continuous space for molecular diffusion limited by the surface boundary. With the aim to simulate long-range intermittent dynamics, only the greatest connected part Ω c of Ω is considered (in the present study, for all samples c t > 0, Ω c percolates through the periodic minimal image). Intermittent Brownian motion was then simulated using the following advanced random walk approach. An interfacial volume is defined as ∂Ω c × x 0 where x 0 = 0.2 pm is an infinitesimal thickness. Diffusion in the pore cavities is described using regular random walk simulations with a bulk-like self-diffusivity D p s estimated from molecular dynamics. When a molecule center of mass reaches ∂Ω c × x 0 , it remains stopped for a time t S distributed according to an exponential probability density function having a first moment t A . After t S , the center of mass is placed at the distance x 0 from ∂Ω c for a new relocation step. The procedure above leads to intermittent Brownian motion where the residence and relocation steps are distributed according to two underlying probability density functions ψ A (t) and ψ B (t) (having t A and t B as first moments). On the one hand, as illustrated in Fig. 3a, the residence times obey a statistics given by: where ψ A (t)dt is the probability that the residence lasts a time between t and t + dt. While the exponential decay in Eq. (3) provides a generic description of the residence time distribution, power-law distributions can be observed in other specific situations such as in media with surface heterogeneity or complex surface dynamics. However, as will be illustrated below, among possible behaviors, the exponential decay is important as it corresponds to a well-defined underlying thermodynamic picture where desorption corresponds to an activated mechanism. Moreover, considering the mapping between microscopic and mesoscopic tortuosities proposed in what follows, it only relies on the mean residence time and not the exact time distribution. On the other hand, ψ B (t) is the bridge statistics which describes the time distribution between a desorption event and the next first re-encounter within the Only data obtained for ε/k B T = 0.5 are shown here for the sake of clarity, but data for other surface/fluid interaction strengths show similar behavior. The data points are the projection of the molecular dynamics tortuosities τ MD onto the Random walk tortuosities τ RW to obtain the residence (t A ) and relocation (t B ) times; i.e. τ MD $ τ 0 RW ½1 þ t A =t B . The schematics on the right illustrate this mapping. In MD simulations, the self-diffusivity is probed by measuring the mean square displacement of the fluid molecules (here, a single molecule is marked in red but the self-diffusivity is averaged over each individual molecule). In RW calculations, the self-diffusion of particles is probed using a specific algorithm which physically accounts for the adsorption and relocation times t A and t B as illustrated in the corresponding schematic. The vertical error bars are the same as in Fig. 2 while the horizontal error bars are their projection onto the dashed lines.
proximal zone ∂Ω c × x 0 . Such generic bridge statistics in confinement is illustrated in Fig. 3b which shows ψ B (t) for methane confined in the disordered sample CS 1.0 with different fluid/surface interactions. On the one hand, after a plateau in the very short time range (≲ps), ψ B (t) decays as a power law ψ B (t)~t −3/2 . On the other hand, in the long time regime, ψ B ðtÞ / expðÀt=t c Þ as strong confinement in the sample cavities introduces a time cutoff t c in the relocation process since every confined molecule eventually returns to the surface within a finite time. This generic behavior for such a finite i.e. confining medium can be described as 40 : where ψ 1 B ðtÞ corresponds to the bridge statistics for a semi-infinite medium (denoted by the symbol ∞). As shown in Supplementary Notes, ψ 1 B ðtÞ can be determined by considering the trajectory of a molecule starting at a distance x 0 from the adsorbing region located in x = 0 and crossing this interface for the first time at a time t 41 : In this equation, the second equality corresponds to the solution in the limit t ) x 2 0 =D p s . Such expressions are valid for a semi-infinite medium where the probability to return to the surface becomes vanishingly small in the long time limit. Figure 3c shows the tortuosity τ RW as a function of the residence time t A as obtained using random walk simulations for the different CS x samples (only data for ε/k B T = 0.5 are shown here for the sake of clarity). The dashed lines in Fig. 3c are RW results obtained by varying t A in a quasi-continuous manner. As expected, the tortuosity can be rescaled as: where τ 0 RW and t B only depend on the specific CS x sample considered. While t B is the typical relocation time, τ 0 RW corresponds to the geometrical tortuosity obtained for a vanishing residence time (t A → 0). As shown in Fig. 3c, projecting the τ MD values obtained by MD (points) onto the ones obtained by RW (lines), i.e. τ MD = τ RW , allows mapping the molecular and mesoscopic tortuosities. This provides a means to estimate for each sample CS x the residence (t A ) and relocation (t B ) times as a function of the fluid/surface interaction strength ε/k B T (values that cannot be assessed using MD for such complex disordered materials). In more detail, t A and t B are such that Considering that τ 0 RW and t B for a given sample and ε/k B T are uniquely defined from the slope and intercept in Fig. 3c, there is only one set (t A , t B ) that verifies τ MD = τ RW . As shown in our previous work 14 , it should be emphasized that t A and t B can be directly estimated from molecular dynamics when simple pore geometries are considered. However, such calculations turn out to be extremely challenging for disordered porous media because the surface/volume decomposition is a complex ill-defined problem. Energy-based criteria such as surface-fluid energy cutoffs or geometrical criteria such as positions to the interface can be used but they rely on arbitrary choices. In contrast, the approach proposed in the present work provides a means to split the complex diffusivity behavior into residence and relocation steps without having to rely on these arbitrary choices.
Bridging molecular/mesoscopic dynamics in disordered media. The residence and relocation times are upscaled parameters which provide a mean to quantitatively bridge the microscopic and mesoscopic dynamics in porous media through the intermittent Brownian motion formalism. Yet, beyond simple mapping procedures like matching molecular and coarse-grained tortuosities, there is a need to establish robust and quantitative physical behaviors for t A and t B . With this aim, the effect of mean pore size d and fluid/surface interaction strength ε/k B T on t A and t B is shown in Fig. 4. In what follows, we first report a molecular model for the residence time t A and then discuss the behavior of the relocation time t B using the formalism of first passage processes.
Residence time t A . Figure 4a suggests that t A follows an activation law for all samples: In this transition state theory, ΔF * corresponds to the free energy barrier that must be overcome by a fluid molecule to escape from the interaction field generated by the host surface. As for 1=t 0 A , it corresponds to the frequency with which the molecule attempts to escape the free energy minimum where it is located. While the activated behavior observed for t A might appear as a surprising result, it can be rationalized through simple thermodynamic arguments. Let us consider a thermodynamic model where the molecule is either adsorbed in the vicinity of the pore surface or in the pore center. As a first-order approximation, it can be assumed that the free energy difference ΔF~N s ε where N s is the number of surface atoms interacting with the fluid molecule. In other words, with this assumption, the free energy of an adsorbed molecule corresponds to the sum of the interaction energies with each neighboring surface atom while the entropy and fluid-fluid interaction contributions are treated as constant. Considering that ΔF * = δΔF with δ ≳ 1 (since the free energy barrier is necessarily larger than or equal to the free energy difference between the adsorbed/non-adsorbed physical states), the scaling in Fig. 4a indicates that α = N s δ. As shown in Fig. 4b, for all samples (i.e. regardless of pore size d), α~3.6 which leads to N s ≲ 3.6. Such a value, which is independent of the considered structure, seems realistic as this corresponds to an underlying molecular picture where an adsorbed molecule interacts with N s~3 to 4 structure atoms. To validate this interpretation, we calculated for all interaction strengths ε/k B T and porous materials CS x , the radial distribution function g(r) between host carbon atoms and methane molecules. The number of local carbon neighbors N c contributing to the free energy barrier involved in the escape time from surface residence was then estimated by integrating g(r) up to the location corresponding to the Lennard-Jones potential minimum r min ¼ 2 1=6 σ, i.e. N c ¼ R r min 0 4πr 2 gðrÞρdr. Considering all structures and interactions strengths, we found <N c > ¼ 3:6 ± 1 which is consistent with the value obtained for α in Fig. 4b. Figure 4c shows that the prefactor t 0 A , which corresponds to the characteristic timescale for activated molecular desorption from the surface, is of the order of~1 ps-a classical value used in transition state theories and nucleation models in dense liquid states. More importantly, t 0 A $ t 0 A;1 ½1 þ γ expðÀd=ξÞ where t 0 A;1 $ 0:07 ps corresponds to the value for infinitely large pores (vanishing confinement). The typical decay length ξ~2.1 Å is of the order of the molecule size σ, therefore indicating that the correction to the escape attempt time is related to the pore size d. This can be understood by the fact that, for a given free energy barrier ΔF * , strong confinement leads to increasing residence times due to the decrease in molecular paths leading to desorption.
Relocation time t B . Figure 4d shows that the relocation time t B scales as t B $ d=D p s . This result is not completely intuitive as it departs from a straightforward estimate obtained using the pore diffusivity D p s and diffusion domain~d, i.e. t B $ d 2 =D p s . Yet, as described quantitatively in what follows, the scaling t B~d can be rationalized by accounting for the fact that the diffusion, i.e. relocation, time within the confining cavities has necessarily an upper bound (due to the finite pore size, each molecule eventually readsorbs to the surface). This constraint, which is at the root of the scaling t B~d , can be quantitatively predicted by introducing a time cutoff t c in the relocation first passage probability ψ B (t). In addition to t c , we also introduce a short time cutoff t 0 as ψ B (t) is necessarily equal to zero for times shorter than the time t 0 needed for a molecule to travel the minimum bridge of extension x min .
As derived in the Supplementary Notes, with the lower/upper time limits t 0 and t c , ψ B (t) simply writes is obtained by writing the normalization condition R 1 0 ψ B ðtÞdt ¼ 1. The first passage distribution for relocation ψ B (t) allows estimating the mean relocation time t B as: As shown in Supplementary Notes, upon inserting ψ B ðtÞ $ exp Àt=t c ð Þ t À3=2 for t > t 0 (0 otherwise) into Eq. (7), it can be shown that: , this expression simplifies as: t c is associated with a geometrical cut-off length r c which indicates the maximal extension of a bridge. r c is of the order of the pore diameter d and can be written as r c = βd, where β~1 is related to the accessible in-pore horizon. Assuming Fickian diffusion upon relocation, we can write t 0 $ x 2 min =2D p s and t c $ β 2 d 2 =2D p s . As shown in Fig. 4d, by assuming that x min is independent of the pore structure, Eq. (8) provides a reasonable description of the observed scaling t B~d with a negative intercept in d = 0. Yet, as detailed in Supplementary Notes, x min can be estimated from the probability density function of the bridge displacement θ(r) where r the is the end-to-end Euclidean distance of a Brownian bridge 42 [see Supplementary Fig. 10a]. With this refined analysis, as shown in Supplementary Fig. 10b, x min does depend on the pore diameter d. Taking into account this dependence, the simulated data t B D p s in Fig. 4b as a function of d can be retrieved using a unique value β0 .7 for all values ε/k B T, as shown in Supplementary Fig. 10c.

Discussion
The statistical physics approach reported in this paper provides an efficient mean to upscale microscopic dynamics in complex porous media to the engineering, i.e., continuum, level. This general and versatile method consists of upscaling molecular constantstypically, the adsorption strength and self-diffusivity-as obtained using molecular dynamics through the formalism of intermittent Brownian motion. While this robust framework is well-established for ordered materials with regular pore geometry and simple pore network topology, the present work extends its scope to ultraconfining disordered porous media with underlying complex freeenergy landscapes. In particular, despite the complex interfacial dynamics in media involving ill-defined surface/volume regions, mapping of molecular dynamics simulations onto intermittent random walk provides a simple yet robust description through the mean surface residence (t A ) and in pore relocation (t B ) times. More importantly, using disordered porous materials with different porosities ϕ/pore sizes d but also fluid/surface interaction strengths ε, t A and t B are found to derive from basic physical models with parameters available to simple experiments. On the one hand, the mean residence time t A is simply related to the fluid/surface interaction strength ε as it corresponds to the characteristic molecular escape time from a low (molecule in the surface vicinity) to a higher (bulk-like molecule in the pore center) free energy state separated by a free energy barrier ΔF *~ε /k B T.
On the other hand, t B can be simply predicted from the confined in-pore self-diffusivity D p s and the corresponding mean-first passage probability distribution which is truncated to account for the finite relocation time in confining cavities. Considering the mesoscopic, i.e., coarse-grained, description adopted in this approach, it is remarkable that all the problem complexity is embedded into two characteristic timescales that are related using simple physical laws to intrinsic material/fluid descriptors. Such upscaling strategy could prove to be useful in numerous fields involving fluid adsorption and transport in porous materials: chemistry (e.g., adsorption, catalysis), chemical engineering (e.g., separation, chromatography), geosciences (e.g., pollutant transport), etc. In particular, among important examples relevant to such practical fields, the present approach can help describe molecular diffusion in the following applications: phase separation of gaseous or liquid effluents through porous media, filtration of small micro-pollutants such as organic/biomolecules, metallic and ionic complexes in water remediation, kinetics of products, reactants and by-products in catalytic processes, etc. From a practical viewpoint, conducting the exact upscaling strategy as reported in this paper can be quite involved; it requires building realistic porous material models and conducting both atom-scale and mesoscopic random walk simulations. However, the physical behavior of t A and t B as established above provides simple rules to predict the long-time fluid diffusion within a given porous material. In practice, all parameters needed to predict this macroscopic behavior are easily accessible experimentally; this includes the pore size d, the fluid/surface energy ε, and the selfdiffusivity D p s . While d can be estimated using adsorption-based techniques or derived using structural data, the fluid/surface energy can be probed from calorimetry or simply estimated from data for similar fluid/solid couples. As for D p s , a good approximation is to take this parameter equal to its bulk counterpart but more accurate data can be estimated by measuring the confined diffusivity using neutron scattering or NMR relaxometry. Inversely, starting from experimentally measured self-diffusivity in confinement, t A and t B can be extracted to shed light on physical phenomena occurring upon fluid adsorption, catalysis, etc. in a given porous material. In this context, our strategy can be coupled with free energy landscape computation to estimate the residence and relocation times. Such calculations are suitable for regular porous materials such as zeolites or metal-organic frameworks (for which dealing with a small porous subspace is sufficient thanks to symmetry considerations). However, such free energy approaches are nearly impossible for disordered porous materials with large representative elementary volume so that an effective approach based on simple physical laws is sound and robust.
Beyond regular adsorption/diffusion processes, our upscaling approach can be used to predict long-time effective diffusivity in problems involving more complex phenomena as observed in natural or anthropic disordered materials (wood, cement, etc.). This includes fluid/solid systems in which desorption is an activated process 43 but also processes involving reactive transport 44,45 and poromechanical effects such as adsorptioninduced swelling 46 . Finally, the present approach can be used to obtain the elementary bricks to be implemented in mesoscopic numerical techniques such as finite elements calculations, pore network models 47 , Lattice Boltzmann simulations but also more formal statistical physics approaches 20,48-50 . As already stated, our mapping procedure is expected to apply to a broad class of fluid/solid couples but some possible limitations must be considered as they can lead to more complex behaviors. Such limitations include the possible role of rich molecular interactions that are potentially long-ranged (e.g. electrostatic). Complex host solids with long-range pore correlations (fractal, low dimension) can also lead to additional complexity. In particular, in extremely narrow pores, confinement induces specific mechanisms such as molecular sieving 51 or single file diffusion 18 that depart from the Fickian regime considered here. Moreover, by considering only percolating matrices (c t > 0), the present study does not address connectivity aspects which can lead to anomalous temperature behavior depending on the ratio of adsorption and connectivity effects 51 .

Methods
Porous material models. Different samples of densities ranging from 0.5 g/cm 3 up to 1.4 g/cm 3 were produced using the following method. For a given density ρ s , the atoms are placed randomly in a cubic box of a size 100 Å (an H/C atomic ratio0 .091 was selected as it corresponds to a typical, realistic value for such disordered porous carbons 52,53 ). Starting from a large temperature, each molecular structure was quenched using molecular dynamics performed using the large-scale atomic/ molecular massively parallel simulator (LAMMPS 29 ). Molecular interactions were described using the reactive empirical bond order (REBO) potential 28 to allow for bond formation/breaking. The quenching procedure is performed in the NVT ensemble by continuously decreasing the temperature from 3000 K down to 300 K in the course of a 5 ns simulation run. Three representative structures are presented in Supplementary Fig. 1 and all .xyz structure files are available upon request.
Grand canonical Monte Carlo. We simulated methane adsorption isotherms at 111.7 K in the various host structures (Supplementary Fig. 2) using Grand Canonical Monte Carlo (GCMC) with the Lennard-Jones parameters gathered in Supplementary Table 2. The saturating vapor pressure of methane at this temperature is P 0 = 101325 Pa (boiling point). In GCMC simulations, we consider a system at constant volume V (the host porous solid) in equilibrium with an infinite reservoir of molecules (methane) imposing its chemical potential μ and temperature T. For a given set T; μ ð Þ, the adsorbed amount is given by the ensemble average of the number of adsorbed molecules versus the pressure P of the gas reservoir (the latter is obtained from the chemical potential according to the equation of state for the bulk gas). The adsorption isotherm is simulated by increasing or decreasing the chemical potential of the reservoir.
The skeleton is considered rigid and the energy U αβ (i, j) between the site i of type α and the site j of type β is given by 54 :  Table 2 for the interactions between sites of the same type, the cross interactions being computed from the Lorentz-Berthelot rules: Molecular dynamics. The methane-saturated structures obtained by GCMC are then used as starting structures for molecular dynamics (MD) simulations. All MD simulations are performed with LAMMPS 29 using the lj/cut potential with the same same-site parameters as the ones used for the GCMC simulations. In all simulations, the porous solid is kept frozen while the probe molecules are simulated at a temperature of 450 K for a NVE production run of 20 ns after a NVT thermalization run of 500 ps. The integration time step is 1 fs and the configurations are saved every 1 ps. To assess the influence of fluid/surface interaction on the effective diffusivity-and, hence, the tortuosity-the fluid/surface interaction strength ε was varied. In so doing, the repulsive interaction felt by the fluid molecules decreases upon decreasing ε so that the porosity explored by the confined molecules increases (inset Supplementary Fig. 3). Consequently, due to this effect, the tortuosity for a given structure strongly evolves with ε without being per se an effect of the interaction strength. To correct this effect, we developed a modified Lennard-Jones potential that keeps the repulsive contribution constant. This modified interaction potential uses a smooth sigmoid function U(r) defined as: UðrÞ ¼ L À ðrÞe sr c þ L þ ðrÞe sr e sr c þ e sr ð11Þ where s = 50 and r c = 0.97σ are the slope and center of the sigmoid, respectively. L − (r) and L + (r) are the connected functions defined for r < r c and r > r c . To keep the repulsive interaction constant, L(r) was maintained fixed as: with ε − = k B T. As shown in Supplementary Fig. 3, upon varying ε + , a modified Lennard-Jones potential with different fluid/surface interaction strengths can be defined while keeping the repulsive part constant.
Topology characterization and diffusion pore space. The pore space Ω available for the dynamics of the spherical methane molecule inside the carbon matrix was determined as follows. A 3D lattice network is first defined with a voxel size 0.02 NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-21252-x ARTICLE nm. A voxel belongs to Ω if its distance to any carbon centers is above 3.605 Å (this value is used as it corresponds to the Lennard-Jones parameter σ for the fluid/ surface interaction). A voxel belonging to Ω is set to 1 (0 otherwise). Such 3D lattice network allows defining the surface boundary ∂Ω of Ω made of surface voxels at the border between Ω and its complementary space. This allows us to define a continuous space for molecular diffusion which is limited by the surface boundary. An interfacial volume is defined as ∂Ω c × x 0 where x 0 is a thickness equal to 0.2 pm. Supplementary Fig. 5 illustrates this procedure by showing for the sample CS 1.0 the resulting digitized pore network and the corresponding retraction graph obtained with a porosity ϕ = 0.177. The molecular trajectory can be described as an alternate succession of a surface adsorption step on ∂Ω c × x 0 followed by a Brownian motion in the confined bulk Ω c leading to a new relocation on the surface. The time step for the Brownian motion is set to 0.1 ps and the selfdiffusion coefficient is estimated from molecular trajectories (mean square displacements) as obtained from molecular dynamics at very early time steps.

Data availability
The data sets and molecular configurations generated during and/or analyzed during the current study are available from the corresponding authors upon request. All MD simulations were performed using the software LAMMPS (stable release from August 31st, 2018