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³—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²⁸ in LAMMPS²⁹ 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 Å.

Fig. 1

Prototypical disordered nanoporous structures.

a Molecular configuration for methane in a disordered nanoporous carbon (here CS_0.70, box length 100 Å). The cyan spheres are methane molecules while the yellow surface represents the carbon porous network. b Mean pore diameter d (top), porosity ϕ (middle), and connection number c_t (bottom) as a function of structure density ρ_s of the digitized pore networks. The same color code applies throughout the manuscript. The samples obtained at larger densities—shown as white/open circles—possess unconnected porous networks (c_t < 0) so that no diffusion is observed in the long time limit.

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–34:

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 self-diffusivity 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)¹⁶. As a result, due to the formal equivalence between permeance and collective diffusivity, i.e. K = D_c/ρk_BT ~ D_s/ρk_BT, the self-diffusivity also provides key insights into transport mechanisms under flow conditions as induced by pressure/chemical potential gradients. The LJ fluid parameters (σ₀ = 3.81 Å, ε₀/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_BT = 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ε₀/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 $$ <\mid \mathbf{r}\left(t\right)-\mathbf{r}\left(0\right){\mid }^2>$$ as a function of time t for methane confined in the different disordered nanoporous materials (only data for ε/k_BT = 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 long-range 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_s – which is inferred from the Fickian regime in the long time limit $$ D_s={\lim }_{t\to \infty }<\mid \mathbf{r}\left(t\right)-\mathbf{r}\left(0\right){\mid }^2>/6t$$ – is lower than the bulk self-diffusivity $$ D_s^0$$. As a result, the tortuosity $$ {\tau }_{\mathrm{MD}}=D_s^0/D_s$$ – defined 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_BT, 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~{\rho }_s^{-x}$$ 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²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 $$ <\mid \mathbf{r}\left(t\right)-\mathbf{r}\left(0\right){\mid }^2>$$ displays a smooth behavior $$ \sqrt{<\mathrm{\Delta }r{\left(t\right)}^2>}~\sqrt{t}$$ from which a confined self-diffusivity can be derived. Yet, Fig. 2b reveals that 4πr²G_s(r, t) displays a complex behavior characteristic of stop-and-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 σ₀, 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³⁵ 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³⁶.

Fig. 2

Complex stop-and-go fluid dynamics.

a Tortuosity $$ {\tau }_{\mathrm{MD}}=D_s^0/D_s$$ as a function of the average fluid/surface energy $$ <U_{fw}>/k_{\mathrm{B}}T$$ as assessed from MD for methane in a disordered nanoporous carbon (same color code as Fig. 1). The dashed lines are linear fits (note that the plot shows $$ \ln \left({\tau }_{\mathrm{MD}}\right)$$). The self-diffusion coefficients D_s are obtained from the mean square displacements shown in Supplementary Fig. 6. The inset shows the in-pore diffusivity $$ D_s^p$$ as a function of mean pore size d for various fluid/surface interaction strengths ε/k_BT. The dashed lines are fits against Eq. (2). The different parameters of this effective fit are shown in Supplementary Fig. 7. The error bars correspond to the standard deviations around $$ {\lim }_{t\to \infty }<\mid \mathbf{r}\left(t\right)-\mathbf{r}\left(0\right){\mid }^2>/6t$$. b Contour plot showing the self-correlation function 4πr² × G_s(r, t) as a function of r and t for methane in a disordered nanoporous carbon (here CS_0.70 with ε/k_BT = 0.5). The local color at a given r, t position indicates 4πr² × G_s(r, t) with the color scale shown on the right. The black dashed lines denote the mean pore size d and mean displacement averaged over all molecules $$ \sqrt{<\mathrm{\Delta }r{\left(t\right)}^2>}~\sqrt{<\mid \mathbf{r}\left(t\right)-\mathbf{r}\left(0\right){\mid }^2>}$$. 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.

To shed more light into the complex diffusivity landscape in such disordered porous media, a single trajectory $$ R\left(t\right)=\sqrt{{\left(\mathbf{r}\left(t\right)-\mathbf{r}\left(0\right)\right)}^2}$$ 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 σ₀. 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 with increasing the fluid/surface interaction strength ε (Supplementary Fig. 8). The second narrow stripe centered at about r ~ 4 Å( ~ σ₀) 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 size – r > d – as 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_s^p$$ can be assessed from the linear regime observed in the short time scale where $$ <\mid \mathbf{r}\left(t_d\right)-\mathbf{r}\left(0\right){\mid }^2>\le 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²/6t_d. The simplest effective framework consists of writing the effective pore-scale diffusivity $$ D_s^p$$ as an average over the whole pore volume, $$ D_s^p=1/N\times \int \rho \left(\mathbf{r}\right)D_s\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}$$ where ρ(r) and D_s(r) are the local density and self-diffusivity at a position r. Within the transition state theory, the bulk self-diffusion coefficient can be written as an activated process $$ D_s^0~\exp \left[-\mathrm{\Delta }{F}^0/k_{\mathrm{B}}T\right]$$ where ΔF⁰ 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_s^0$$ 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⁰ − ζ(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_{\mathrm{s}}\left(r\right)=D_{\mathrm{s}}^0\exp \left[\zeta \left(r\right)/k_{\mathrm{B}}T\right]$$¹⁹. 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: $$ D_s\left(r\right)=D_s^s$$ for r > d/2 − σ while $$ D_s\left(r\right)=D_s^0+\left(D_s^s-D_s^0\right)\exp \left[-\left(d/2-\sigma -r\right)/r_0\right]$$ for r ≤ d/2 − σ (where $$ D_s^s$$ is the surface self-diffusivity in the vicinity of the pore surface while $$ D_s^0$$ 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₀ towards the bulk diffusivity $$ D_s^0$$ 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:

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³⁹. Moreover, such surface/volume partition and the resulting predictions in terms of in-pore diffusivities $$ D_s^p$$ 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₀ where x₀ = 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_s^p$$ estimated from molecular dynamics. When a molecule center of mass reaches ∂Ω_c × x₀, 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₀ 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:

Fig. 3

Mapping microscopic dynamics onto mesoscopic intermittent brownian motion.

(a) Typical residence statistics ψ_A(t) with t_A = 0.1 ns. The dashed line shows the expected scaling $$ {\psi }_A\left(t\right)\propto e^{-t/t_A}$$. (b) Bridge statistics ψ_B(t) for CS_1.0 with two fluid/surface interactions ε/k_BT as indicated in the graph. The solid lines are fits against Eq. (4), with t_c = 42 ± 0.5 ps for ε/k_BT = 0.3 and t_c = 51 ± 0.7 ps for ε/k_BT = 0.5. The dashed lines indicate the main regimes: ψ_B(t) ∝ t^−3/2 and $$ {\psi }_B\left(t\right)\propto e^{-t/t_c}$$. (c) Random walk tortuosity τ_RW as a function of the residence time t_A (dashed lines). Only data obtained for ε/k_BT = 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. $$ {\tau }_{\mathrm{MD}}~{\tau }_{\mathrm{RW}}^0\left[1+t_A/t_B\right]$$. 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.

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_BT = 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:

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_BT 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.

Fig. 4

Residence and relocation times.

a Logarithm of the residence time t_A as a function of ε/k_BT where t_A is obtained by mapping the molecular and mesoscopic tortuosities as shown in Fig. 3c. The dashed lines are fits in the form $$ t_A=t_A^0\exp \left(\alpha \epsilon /k_{\mathrm{B}}T\right)$$ with $$ t_A^0$$ the characteristic residence time for vanishing fluid/surface interactions. The error bars are those given in Fig. 3c. b Impact of pore size d on the dimensionless parameter α; α = 3.6 ± 0.4 ∀ p. The error bars correspond to fitting errors. c Characteristic residence time $$ t_A^0$$ as a function of mean pore size d. The dashed line is a fit against an effective decaying function $$ t_A^0\propto t_{A,\infty }^0\left[1+\gamma \exp \left(-d/\xi \right)\right]$$ with the critical pore size ξ = 2.1 ± 0.5 Å and $$ t_{A,\infty }^0~0.07$$ ps. d Relocation time $$ t_B\times D_s^p$$ as a function of the mean pore size d for ε/k_BT = 1. As expected, $$ t_B\propto d/D_s^p$$ as the finite time spent upon relocation in the cavities introduces an upper bound in the survival probability (see text).

Relocation time t_B

Figure 4d shows that the relocation time t_B scales as $$ t_B~d/D_s^p$$. This result is not completely intuitive as it departs from a straightforward estimate obtained using the pore diffusivity $$ D_s^p$$ and diffusion domain ~ d, i.e. $$ t_B~d^2/D_s^p$$. 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₀ as ψ_B(t) is necessarily equal to zero for times shorter than the time t₀ 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₀ and t_c, ψ_B(t) simply writes $$ {\psi }_B\left(t\right)=C\times \exp \left(-t/t_c\right)/t^{3/2}$$ for t ∈ [t₀, t_c] where $$ C={\left[2\exp \left(-t_0/t_c\right)/\sqrt{t_0}-2\sqrt{\pi /t_c}\mathrm{erfc}\left(\sqrt{t_0/t_c}\right)\right]}^{-1}$$ is obtained by writing the normalization condition $$ {\int }_0^{\infty }{\psi }_B\left(t\right)\mathrm{d}t=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 $$ {\psi }_B\left(t\right)~\exp \left(-t/t_c\right)\times t^{-3/2}$$ for t > t₀ (0 otherwise) into Eq. (7), it can be shown that: $$ t_B=C\times \sqrt{\pi t_c}\mathrm{erfc}\left(\sqrt{t_0/t_c}\right)$$. By writing that t₀ ≪ t_c (i.e. $$ C~\sqrt{t_0}/2$$), 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_{\min }^2/2D_s^p$$ and $$ t_c~{\beta }^2{d}^2/2D_s^p$$. 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⁴² [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\times D_s^p$$ in Fig. 4b as a function of d can be retrieved using a unique value β ~ 0.7 for all values ε/k_BT, as shown in Supplementary Fig. 10c.

Porous material models

Different samples of densities ranging from 0.5 g/cm³ up to 1.4 g/cm³ 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 ratio ~ 0.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²⁹). Molecular interactions were described using the reactive empirical bond order (REBO) potential²⁸ 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₀ = 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 $$ \left(T,\mu \right)$$, 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⁵⁴:

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²⁹ 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:

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 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₀ where x₀ 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₀ 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 self-diffusion coefficient is estimated from molecular trajectories (mean square displacements) as obtained from molecular dynamics at very early time steps.