Tuning the permeability of regular polymeric networks by the cross-link ratio

Nano-, micro-, and mesoporous materials are promising candidates for numerous applications, such as particle absorption, storage, drug delivery, molecular sieving, and catalysis.^1–11 Nowadays, there exists a big variety of composites in engineering, making up a range from soft to solid materials: Many of them contain metals, while others are purely organic. Several natural microporous materials, such as zeolites, are based on hybrid metallic–organic composites, and these own a rather random microstructure. They have inspired the design and synthesis of the so-called metal–organic frameworks (MOFs) and porous coordination polymers¹² (PCPs) featuring periodic structures with well-defined network mesh sizes. In terms of porous organic chemicals, synthetic polymer materials have undergone a similar evolution^13–15 and provide a wide range of polymer networks from rather soft and randomly cross-linked ones, such as hydrogels^16–19 or aerogels, to highly ordered and much more rigid covalent organic frameworks^20–25 (COFs), porous aromatic frameworks^26,27 (PAFs), and (micro-) porous polymer networks^28,29 (PPNs) or high-strength hydrogels, such as tetra-arm polyethylene-based networks.^30–34 In fact, this list is not complete and the intense ongoing research in polymeric materials raises the expectation that arbitrary and regular network architectures with desired mesh sizes and properties can be fabricated.

Although the physical and chemical properties of all these materials do differ, they are essentially networks of polymers with interconnecting junctions. We will denote these junctions as cross-linkers throughout this work. By identifying and counting the chain-monomers and the linkers, all networks can be characterized by the cross-link ratio expressed as

with N_x and N_mer being the cross-linker number and chain-monomer number in the polymer network, respectively. For convenience, we denote the chain-monomers simply as monomers in the following.

The polymeric networks are, in adaptation to their function, manufactured at various cross-linker concentrations, chain lengths, mesh sizes, or volume fractions. Regarding stimuli-responsive hydrogels, the cross-link ratio is usually below 15 percent in order to retain its responsiveness.³⁵ However, highly cross-linked gels and many frameworks (MOFs, COFs, and PAFs) used as absorbing and storage devices show values of α ranging from 15% to 50%. The cross-links reinforce the network rigidity^31,36 and, in addition, affect the uptake and the mobility of solutes in the network,^1,37,38 which is the focus of the present work. In fact, the permeability of the network to solutes, which we interchangeably call tracers or penetrants, is of key importance to smart material design^39–44 since its control allows the fine-tuning of the function and performance of devices such as drug carriers, filtration membranes, and nanoreactors.

In many experiments and computer simulations, a planar slab-like geometry of the network is very convenient since the analysis can be carried out in one dimension. The permeability across such membranes can be obtained by the solution–diffusion model, which is based on Fick’s first law of diffusion.^45,46 With the assumption of very small concentrations and their gradients, the linear response limit is applicable, and the steady-state solutions equal the equilibrium situations of such membranes. The permeability can thus be calculated by the product of the equilibrium partition ratio K (a measure for the uptake) and the tracer self-diffusion inside the network D_in (quantifying the mobility), reading^{39,40,43,47–52}

Note that the definition of the permeability may alter with respect to the subject of research and the definition of the partition (or sorption, or solubility) coefficient/ratio. We have decided to regard P as a material property that depends on the type of membrane and the solutes, i.e., it is independent of the membrane width, and has the units of the diffusion coefficient. We know from our group’s previous work^53,54 that the permeability exhibits non-trivial functions of the polymer volume fraction and the interaction strengths between the tracers and the network. These relations may be non-monotonic as a result of the possibly anti-correlated behavior of the tracer partitioning and diffusivity. Maximizing or minimizing the solute permeability can be tuned by the polymer volume fraction, which is controlled by the cross-link ratio α, the key parameter in the present work.

In this work, we emphasize another important influence on the permeability, namely, a dissimilar chemistry of cross-linkers and monomers. In combination with the solute and solvent chemistry, various complex local interaction scenarios arise, which are caused by hydrogen bond formations, hydrophobic (or other entropic) effects, and (local) electrostatic and steric contributions.⁵⁵ By quantifying the solute adsorption to specific parts of the polymer, i.e., the chain and the cross-linkers, the complexity on the atomistic level can be reduced and mapped onto coarse-grained potentials with effective interaction parameters. This enables us to systematically vary the effective interaction strengths between the tracers and the cross-linkers as well as between the tracers and the monomers.

We expect various non-monotonic relations between the cross-link ratio and the partitioning and consequently an even more complex permeability. Similar non-monotonic behaviors of the partitioning and diffusivity in dependence on the chemistry of hydrogel and the solutes have been observed experimentally for different degrees of co-polymerization and for different solute types.^56,57 Although many numerical studies investigate the influence of cross-linkers, polymer density, many-body effects, particle size, and qualitatively different interactions, on the solute partitioning and diffusion,^53,58–71 none of them, to the best of our knowledge, has introduced an independent, specific interaction strength between cross-linkers and the solutes, as proposed herein.

Our investigation shows how the cross-link ratio may be used to tune the membrane’s solute permeability when the polymer network of the membrane retains a regular architecture and conformation. We provide insight into the intricate interplay between partitioning and diffusion in dependence on the tracer–polymer interaction strengths. To this end, we perform coarse-grained, Langevin dynamics simulations (implicit solvent) of polymer network membranes with different cross-link ratios and tracer–polymer interaction strengths. The cross-linkers are tetrafunctional, and the polymer chains are of equal length, and hence, the network owns a tetrahedral (diamond-like) architecture (see Fig. 1 for simulation snapshots and Fig. 2 for an illustration of different regular architectures). We employ Lennard-Jones (LJ) potentials for non-bonded interactions and use harmonic potentials for bonds and angles. For the generality of our model and transferability to various porous materials, we present our results in the standard (LJ) reduced units.

In our coarse-grained simulations, all particles are considered spherical and equal in size and mass. The effects of the solvent, i.e., the numerous collisions of the particles with the solvent molecules, are implicitly imposed by dissipation (friction) and fluctuating forces, as described by the standard Langevin dynamics framework. The thermal fluctuation of the ith particle is modeled by a Gaussian white noise $ξil(t)$ (with direction l ∈ {x, y, z}). It is delta-correlated in time, independent of any other particle j, and orthogonal in space, reading $⟨ξil(t)ξjl′(t′)⟩=2kBTγδ(t−t′)δijδll′$⁠. The noise intensity k_BTγ is proportional to temperature T and the Boltzmann constant k_B and is linked to the friction coefficient γ and related to D₀, the unperturbed diffusion of suspended particles in infinite dilution, via D₀ = k_BTγ⁻¹.

The non-bonded particle interactions are subject to pairwise LJ potentials of the form

with r_ij = |r_i − r_j|, σ being the LJ length (regarded as the particle diameter) at which U_LJ vanishes, and $ε$_ij being the potential depth at 2^1/6σ, which depends on the types of the particle pair ij. In the simulations, the LJ interactions are truncated at the standard cut-off length 2.5σ.

Our system consists of three particle types, namely, “monomers,” “cross-linkers,” and “tracers.” The parameter $ε$_mer corresponds to the tracer–monomer interaction, while $ε$_x denotes the tracer–cross-linker interactions. They will be altered in the range from 0.01k_BT to 3.0k_BT and from 0.01k_BT to 2.0k_BT, respectively. All remaining potential depths (tracer–tracer and polymer–polymer) are fixed at $ε$ = 0.1k_BT (repulsive/good solvent regime⁷³). For convenience, we present the interaction strengths in reduced energy units, i.e., as β$ε$⁠, with $β=(kBT)−1$⁠.The bonded interactions (bonds and angles) in the polymer network are modeled by harmonic springs, reading

The bond potential U_b is identical for all pairs with κ_b = 100k_BT/σ² and r⁰ = σ. Concerning the angle potential, all constants are equal (κ_a = 10k_BT/rad²), but we distinguish between two types of equilibrium angles $ψj0$⁠, determined by the angle’s vertex particle j. The vertex can be either a monomer (⁠$ψmer0=π$⁠) or a cross-linker (⁠$ψx0=arccos(−1/3)=θ≈109.47°$⁠, the tetrahedral angle). The equilibrium angles and relatively high spring constants keep the polymer network close to the regular extended shape. The chosen values for the angle stiffness κ_a and the number of tracers ensure our dilute concentration regime in consideration, which is further discussed in Sec. SI of the supplementary material.

We construct polymer network membranes, which are periodic in two dimensions and with a finite thickness in the z-direction. In order to study the solute partitioning and their diffusion, we set up systems with two distinguishable phases, i.e., the bulk reservoir phase and polymer network phase. The size of the polymer network and the bulk phase are chosen large enough to minimize interfacial effects. The simulation box length in x-, y- and z-directions are denoted as L_x = L_y = L_xy and L_z, respectively. The actual values depend on the cross-link ratio α, discussed below, but have a fixed aspect ratio of L_z = 3L_xy. Simulation snapshots of equilibrated networks with three differing cross-link ratios are depicted in Fig. 1.

The cross-linkers are initially placed on regular diamond lattice sites, and (nearest) neighboring cross-linkers are interconnected by polymer chains of equal lengths, i.e., equal number of monomers N_chain per chain, which is related to the cross-link ratio of a fully periodic network as $α=(2Nchain)−1$⁠. This work studies networks with nine different cross-link ratios, summarized in Table I. The maximum number of monomers between two cross-linkers was chosen to be N_chain = 25, corresponding to a very low cross-link ratio of α = 0.02. The shortest distance between two cross-linkers is N_chain = 1, resulting in a relatively high α = 0.5. Owing to our setup procedure, it is convenient to express the size of the initial polymer network by the amount of diamond lattice unit cells. One membrane has the same number of unit cells in each direction (without accounting for the x–y-periodicity). The least dense network (N_chain = 25, α = 0.02) consists of exactly one unit cell. The number of unit cells of the networks with higher α is chosen such that the width L_xy is similar to that with α = 0.02 (see Table I). After the initial placement of cross-linkers and chain-monomers, the network (without tracers) is equilibrated with the stochastic Langevin integrator with γ/m = τ⁻¹ (⁠$τ=σβm$⁠) together with an isotropic Berendsen barostat⁷⁴ (τ_p = 0.1τ) with zero equilibrium pressure. This NpT equilibration is carried out for 5 × 10³τ (10⁶ time steps with one increment of dt = 0.005τ) to ensure the relaxation of the polymer network and the convergence of the simulation box volume, i.e., the network is initially stress-free.^60,61

In the subsequent step, tracer molecules are added randomly into the simulation box. The number of tracer particles N_tracer ensures a relatively high dilution and sufficient sampling. The average concentration yields $c¯=Ntracer/(Lxy2Lz)≈0.005/σ3$⁠, yet testing with a higher dilution, i.e., $c¯≈0.001/σ3$⁠, has not significantly changed the results (see Sec. SI of the supplementary material).

The network sizes and the corresponding polymer and tracer particle numbers are summarized in Table I. The tracers and the network are further equilibrated in the NVT ensemble for 5 × 10⁴τ (10⁷ time steps). The final production run (NVT) is carried out for 1.5 × 10⁵τ (3 × 10⁷ time steps). This procedure has been performed for all α and possible interaction parameter combinations of β$ε$_x ∈ {0.01, 0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 1, 2, 3} and β$ε$_mer ∈ {0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.56, 0.5, 0.6, 0.7, 1, 2}. The network creation and the simulations are performed with the LAMMPS^75,76 software package.

Every 100 time steps of the production run, the z-positions of all particle types (polymer, cross-linker, and tracer) are collected to create a total of 3 × 10⁵ particle probability density profiles (particle average) for each type. The center of mass (COM) in the z-direction of all polymer beads for each distribution is determined and is used as a reference point for the analysis, as depicted in Fig. 1(c). 100 subsequent profiles are further averaged (temporal binning), reducing the total number of density profiles per particle type to 3 × 10³.

For the estimation of the membrane width L (along z), we use the full width at half maximum (FWHM) of the polymer number density profiles c_p(z). This construction is illustrated by red dashed lines in Fig. 1(c) (note that we instead present the local polymer volume fraction in this figure, which is proportional to the polymer number density). The membrane width L is used as a measure to divide the system into three domains, namely, the “network” phase (|z| ≤ L/4), the “bulk” phase (L < |z| ≤ L_z), and their “interface” (L/4 < |z| ≤ L). Owing to good sampling, the latter is chosen relatively generous to minimize the polymer effects on the bulk phase. In the network region, we estimate the average polymer volume fraction, i.e., the fraction of the polymer volume V_p in a given volume V, yielding

in which the polymer beads are approximated as hard spheres with diameter σ. Furthermore, the tracer profiles are used to determine the partition ratio as

where c_in and c₀ are the tracer concentrations inside the network and the bulk region, respectively. One example profile with a moderate cross-link ratio of α = 0.1 and attractive interaction strengths (β$ε$_mer = β$ε$_x = 1) is presented in Fig. 1(c).

For the calculation of the diffusion from the mean squared displacement (MSD), the trajectories of the tracer particles are required. The positions of the last third of the production run are output every 1000 time steps (5τ), and the z-coordinates are evaluated relative to the z component of the polymer’s COM. The particle trajectories are split into parts (referred to as segments k) of length t_seg = 500τ that are inside the polymer network phase (see Sec. SII A of the supplementary material for more details). With t_seg ≫ 1/γ, we ensure to measure the long-term diffusion for which the MSD ∝ t holds. The segments were used to calculate the tracer particle diffusion in x- and y-directions inside the network, yielding the 3D effective diffusivity as

where ⟨⋯⟩_seg denotes the average over trajectory segments of all particles k identified as inside the network phase and $Dxy2D$ denotes the effective diffusion in the x- and y-directions only. Because of the symmetry of the network, we assume all one-dimensional diffusion coefficients (in x, y, and z directions) to be identical, decoupled, and additive, as it is the case for the free diffusion with orthogonal stochastic forces $ξil(t)$⁠. This assumption has been verified by the 3D diffusion result from simulations of fully periodic networks. A more detailed description of the data preparation, the calculation of Eq. (8), the verification of this method, and examples of the time resolved MSD can be found in Sec. SII A of the supplementary material.

Although the present work focuses on regular networks with a diamond-like architecture, the analytical results for the penetrant partitioning, diffusion, and permeability presented in Secs. III B–III D can be extended to different network topologies if the relationship between the cross-link ratio α and the polymer volume fraction ϕ_p is known. Furthermore, ϕ_p is a key parameter that tunes the partitioning and diffusion and it is affected by the architecture of the network (expressed by the cross-linker functionality f) as well as the cross-link ratio α.

We first consider regular and fixed (no fluctuations) networks in the extended state. More precisely, the cross-linkers are placed on lattice sites and are interconnected by straight chains of equal length (N_chain monomers). All entities, i.e., monomers and cross-linkers, have the same diameter σ and equal distance r⁰ to bonded neighbors. For convenience, we will denote such networks as “static” throughout this work. We evaluate three different regular architectures, which differ in the cross-linker functionality f. These are the simple cubic (f = 6), the body-centered cubic (f = 8), and the diamond lattice (f = 4). The cross-link ratio is hence defined as α = 2/(fN_chain). The reader is referred to Fig. 2 for the illustration of the different architectures and the corresponding polymer volume fractions $ϕpstatic$ as functions of the cross-link ratio α. The mathematical details can be found in Sec. SIII of the supplementary material. We calculate $ϕpstatic$ for all networks without overlapping beads (i.e., distance/bond length equals diameter r⁰ = σ), and for f = 4 we additionally show the effect of overlapping beads (r⁰ = σ/2).

Comparing the different architectures (non-overlapping), $ϕpstatic$ increases with f. However, the scaling with respect to α does not depend on f. For very small α, $ϕpstatic$ scales with α². With values of α just below 0.5, $ϕpstatic(α)$ shows a trend toward linearity. For the overlapping case, we choose r⁰ = σ/2. This case can be considered as the excluded volume fraction $ϕexstatic$ for solute particles with identical diameter as the polymer beads. With f = 4, the scaling with α is similar since their relation is almost linear (⁠$ϕexstatic≈5ϕpstatic$⁠).

Evaluating the flexible model from the simulations (black symbols in Fig. 2), we see that the actual polymer volume fraction ϕ_p is always greater than the static reference $ϕpstatic$⁠. This is clear since the polymer chains are subject to thermal fluctuations and the network is no longer completely extended, decreasing the network mesh/volume. This leads to the fact that α² for low cross-link ratios does not hold. Note that with lower angle stiffness (κ_a) or worse solvent quality (⁠$ε$⁠), ϕ_p will deviate more from $ϕpstatic$⁠.

An additional influence on ϕ_p in the simulation may arise from the tracer–monomer interactions. If $ε$_mer is high and the polymer rigidity is low (i.e., low α) and if the concentration of the tracers (⁠$c¯$⁠) is high enough, a collapse of the polymer network can be induced (see Sec. SI of the supplementary material for a comparison of different concentrations). This “solute-induced” collapse is already described in our previous work.⁶⁹ Although α controls the occurrence of the collapse, ϕ_p will be rather independent of α if the network is in a densely packed configuration. The collapsed cases are not considered in Secs. III B–III D since the analysis of the partitioning and diffusivity is fundamentally different.

In this section, we first describe the simulation results of the partition ratio and subsequently make use of two common analytical approaches that acceptably describe the simulation results. The partition ratio, K_sim = c_in/c₀, shows different qualitative scenarios in regard to the interaction strengths $ε$_x and $ε$_mer and the cross-link ratio α. In Figs. 3(a)–3(c), we show a selection of K as functions of the cross-link ratio representing the tested parameter range.

For equal interaction strengths (⁠$ε$_mer = $ε$_x), K can be monotonically increasing (⁠$ε$_x = $ε$_mer > $ε$^c) or monotonically decreasing (⁠$ε$_x = $ε$_mer < $ε$^c), depending on the critical interaction strength $ε$^c. The critical interaction strength for non-bonded beads in high dilution is known,⁷³ β$ε$^c ≈ 0.32, yet it is expected to deviate slightly for bonded beads because of many-body effects. The cases with equal interaction strengths have been thoroughly studied by our group for comparable polymer membranes^53,69 and the results are qualitatively consistent.

For non-equal interactions (⁠$ε$_mer ≠ $ε$_x), the monotonous scenarios for K(α) are, as one would expect, also observed with $ε$_x, $ε$_mer > $ε$^c or $ε$_x, $ε$_mer < $ε$^c, respectively. By comparing two parameter pairs with interchanged values, e.g., orange symbols in Fig. 3(b) (β$ε$_x = 1.00 and β$ε$_mer = 0.50) vs purple symbols in Fig. 3(c) (β$ε$_x = 0.50 and β$ε$_mer = 1.00), one can recognize that the tracer–monomer interaction strength $ε$_mer has a greater impact on the partitioning since N_mer is always greater than N_x. On some occasions, we observe non-monotonic behavior of K with possible extrema. These can occur when one interaction strength is very low, while the other has moderate to high values. For instance, the orange symbols in Fig. 3(a) (β$ε$_x = 1.00 and β$ε$_mer = 0.01) show a minimum at α = 1/6, and the blue symbols in Fig. 3(b) (β$ε$_x = 0.01, β$ε$_mer = 0.5) maximize at α = 1/4. Since these extrema are hardly visible in the plotted range, see Sec. SIV of the supplementary material for more details. However, these extrema are not very pronounced and the partitioning does not deviate much from unity.

The simulation results are now compared with two theoretical approaches. Both are based on the assumption of an infinite tracer dilution. We first consider static, fully extended networks as discussed above. Hence, the positions r_i of all polymer beads (monomers and cross-linkers) for all α are known, and the potential landscape H(r_j) for a single tracer at any position r_j inside one network unit cell can be computed. This is the superposition of the potentials [Eq. (3)] of all polymer beads in one unit cell and its first periodic images (including the diagonals) within the cut-off (2.5σ), reading

with r_ij = |r_i − r_j| and $ε$_ij depending on i, the polymer type (⁠$ε$_x or $ε$_mer). The probability of finding a tracer at r is proportional to its Boltzmann factor exp(−βH(r)), and the partition ratio is eventually obtained by averaging these probabilities in one unit cell as

The partitioning described by Eq. (10) is in good agreement, yet it overestimates the simulation results [see Figs. 3(a)–3(c)] since neither the polymer conformation nor possible steric self-exclusion (tracer–tracer interactions) is considered by this approach. The testing of different concentrations (see Sec. SI of the supplementary material) shows that the self-exclusion can lower K at most by 10% at $c¯≈0.005σ3$ and only for very high cross-link ratios and interactions strengths $ε$ = 2.0. Furthermore, since the static networks (used for calculating K_μ) have smaller polymer volume fractions than the corresponding flexible networks in the simulation (see Fig. 2), one would expect K_μ to underestimate K_sim. We thus conclude that the differences between K_sim and K_μ stem from non-trivial contributions of the polymer’s local conformation.

Our second theoretical approach applies the virial expansion up to the second order. More precisely, two independent second virial coefficients $B2x$ and $B2mer$⁠, accounting for the two-particle interactions of tracers with cross-linkers and the monomers, respectively, are calculated. These enter the partition ratio as

where the prefactors are derived from the amount of cross-linkers n_x = 8 and monomers n_mer = 8α⁻¹ per unit cell with α-dependent volume V_cell. A linear approximation of this approach was already established in our previous research⁵⁵ by making use of the solute–polymer adsorption coefficients (Γ^* ≈ − 2B₂) obtained from atomistic simulations, which predicted the different qualitative scenarios for the partition ratio of regular networks as a function of α presented herein. The second virial coefficients in the present work are determined by

In the case of the Lennard-Jones interactions [Eq. (3)], the second virial coefficient B₂ vanishes at the critical value β$ε$^c ≈ 0.32, and the partitioning [log(K) ∝ −B₂] thus equals unity.

The partition ratio given by Eq. (11) is compared with the simulation results in Figs. 3(a)–3(c). The B₂ approach overestimates the partition ratio, especially for high interactions strengths $ε$_x, but is in good qualitative agreement. In fact, the monotonicity and even the existence of extrema of K(α) can be predicted. The deviations from the simulation results stem from the many-body effects of the polymers, as these are disregarded by the second virial coefficient.

In this section, we present the simulation results for the tracer diffusion and discuss the limits in which the existing theories are applicable. In Figs. 3(d)–3(f), the normalized penetrant diffusion coefficient D_in/D₀ inside the polymer network vs the cross-link ratio α for a selection of representative parameter sets (⁠$ε$_x, $ε$_mer) is presented. The normalized diffusion is, as expected, always below unity and decreases with the increase in the cross-link ratio at moderate interaction parameters. This can be rationalized by steric obstruction and additional trapping in local potential minima, reducing the probability of free diffusion. Regarding the trapping, the tracer–monomer interaction $ε$_mer has the strongest effect since it contributes to the solute trapping to the chains as well as to the cross-link region. The tracer–cross-link interaction parameter $ε$_x, however, plays a subordinate role. It shows the strongest influence for moderate values of α ∈ [0.1, 0.25].

For strong interactions, e.g., β$ε$_mer = 1 and β$ε$_x = 3, D_in at first decreases and then increases again, exhibiting a minimum at α = 0.25. A more pronounced minimum is found for higher attractions, e.g., β$ε$_x = β$ε$_mer = 2.0, visualized in Fig. S2 of the supplementary material. This effect has already been observed for particle diffusion in static and perfectly periodic media.⁷⁷ Therein, it is supposed that the overlap of potential wells facilitates transport pathways, which speed up the diffusivity. In other words, the overlap flattens the energy landscape, reducing the particle trapping. The diffusion thus tends toward the free-volume limit. Noteworthily, at α = 0.5, the diffusion appears almost independent of $ε$_x.

Now, we briefly overview available theories for the diffusion. In the last decades, many empirical laws and theoretical frameworks have been established that provide some explanation of the dynamics of penetrants in porous media, lattice structures, and flexible polymeric media.^78–82 The concepts scale the diffusion to pore/mesh sizes^83–85 and to free volume and obstruction theories^{41,78,86–89} and may also account for hydrodynamic effects.^89–91 Peppas and Reinhart proposed a theory that accounts for cross-linkers, i.e., they account for the average mesh size in swollen hydrogels.^92,93 In fact, many of these theories describe well the diffusion in the purely repulsive (steric exclusion) limit (⁠$ε$_x = $ε$_mer ≈ 0.01). For instance, a variation of the obstruction theory,⁸⁷ $Din/D0=(1−b1ϕp)2/(1+b1ϕp)2$⁠, the free-volume theory by Yasuda et al.,^41,86

the general hydrodynamic theory,⁹¹ $Din/D0=exp(−b3ϕpb4)$⁠, and even a linear approximation, D_in/D₀ = 1 − b₅ϕ_p, can acceptably be fit (with fitting parameters b_i), where the free volume theory Eq. (13) (with b₂ set to 1 and $bϕYasuda≈5.5$⁠) yields the best description. However, all approaches mentioned so far are only of steric nature. They disregard the attractions in H(r), generated by the numerous polymer beads, whose effects become evident for higher interaction strengths.

Moreover, none of these theories accounts for the possible specific solute–network interactions, which we have introduced with the parameters $ε$_x and $ε$_mer. In previous simulations with very flexible and polydisperse polymeric membranes,⁵⁴ we proposed an extension to the free-volume theory that scales the prefactor b₂ in Eq. (13) with $ε$_p, the tracer–polymer interaction strength. This theory can only be acceptably fit to the simulation data of the present work if α is fixed. We thus provide a new theory that describes the simulation results while accounting more accurately for α and the interaction parameters $ε$_mer and $ε$_x. The interested reader can find the detailed derivation in Sec. SVII of the supplementary material.

The tracer diffusion is considered as a rate (or jump) process within the polymer matrix, and it is governed by the adsorption and desorption at binding sites (traps), well described by Berg for intracellular diffusion in 1986.⁹⁴ It can be formally written as

where τ_free and τ_trap describe the average time the penetrant freely diffuses and the mean time of being trapped, respectively. This result can be interpreted as $Din=D0Wfree$⁠, with $Wfree$ being the probability that the particle freely diffuses, which can be identified in many theoretical approaches.⁸¹ In the steric exclusion limit, i.e., β$ε$_mer = β$ε$_x ≈ 0.01, when no adsorption is at play, $Wfree$ is described by the (free-volume, obstruction, and hydrodynamic) theories mentioned above. For stronger interactions, the trapping of tracer molecules dominates the dynamics, and τ_trap can be expressed in terms of (inverse) Kramers’ rates. The time τ_free to find such a trap has been extensively studied in random and ordered media.^95–98 It is essentially a mean first passage time problem, and its solution depends on the unperturbed diffusivity D₀ and the polymer volume fraction ϕ_p.

The final result can be written in terms of our input parameters α, $ε$_mer, and $ε$_x and the polymer volume fraction ϕ_p(α), reading

where b_ϕ, b_mer, and b_x are dimensionless fitting parameters, which depend on the polymer’s local conformation. The fitting parameters b_mer and b_x scale the contributions (many-body effects) to the specific traps, and b_ϕ is a dimensionless steric prefactor. The binding to the polymer chains only depends on $ε$_mer, whereas for the binding to the cross-link regions, both interaction parameters need to be considered. Be reminded that the polymer volume fraction is a function of the cross-link ratio. It depends on the flexibility of the polymer chains, the network architecture (Fig. 2), and the solvent quality and is thus kept as an empirical input parameter. Noteworthily, with only one trap type (⁠$ε$_x = $ε$_mer), Eq. (15) is formally similar to the diffusion model for molecules in semidilute polymer solutions proposed by Petit et al.⁹⁹ Although the derivation and argumentation are somewhat different, the physical concept, i.e., regarding diffusion as a hopping process, is similar to Berg’s initial idea.

The final analytical model [Eq. (15)] has been fit against 880 out of 1170 simulation results with parameters sets (⁠$ε$_x, $ε$_mer) at once, i.e., we obtain unique results for $bϕ,bmer$⁠, and b_x. As mentioned, only small to moderate values of the cross-link ratio (α ≤ 0.25) were taken into consideration since at α = 0.5, the local trapping is reduced and the assumptions are violated. Furthermore, the most sparse networks with α = 0.02 were also excluded as a consequence of large uncertainties and unacceptable aspect ratios of the membranes (see Table I), violating the symmetry argument for $Din=Dxy2D$⁠. In the case of a network collapse due to strong tracer–monomer interactions, the assumptions of our theory do not hold, and we therefore omitted all simulation results with β$ε$_mer = 2 at $α∈0.025,0.033,0.05$⁠. The inverse errors of the $Dinsim/D0$ have been used as fitting weights. The results of the fit are b_ϕ = 6.6 ± 0.5, b_mer = 0.9 ± 0.1, and b_x = 0.65 ± 0.05. The errors are estimated from a comparison with an unweighted fit of the same data.

The fit is in overall good agreement and allows an interpretation based on the parameter values. The dimensionless steric parameter b_ϕ accounts for the effective surface area per polymer volume fraction and the characteristic length scale of the free diffusion. It is specific to the type of polymer, its network architecture, the solvent, and the tracer molecules. Nevertheless, for very small α (ϕ_p → 0) and zero attractions (⁠$ε$_mer = $ε$_x to 0), one can show that Eq. (13) [or Eq. (S19)] and Eq. (15) converge to D_in/D₀ = 1 − b_ϕϕ_p. Thus, the steric limit described by the free volume theory is incorporated in our theory for small α. The parameters b_mer and b_x describe the contribution of $ε$_mer and $ε$_x to the specific traps. As already observed in the simulation results, the tracer–monomer interaction has a greater impact on the diffusivity than the interactions between tracers and cross-linker beads, even for the traps located at these junctions. Since b_mer ≈ 1, the trapping approach appears to be meaningful, and the value of b_x ≈ 0.65 makes sense since the cross-linker is somewhat shielded from the surrounding monomers, lowering the contribution of $ε$_x to this trap.

The simulation results and two theoretical approaches of the solute permeability (P = D_inK) are depicted in Figs. 3(g)–3(i) as a function of the cross-link ratio. In general, as already observed for the partitioning and diffusion, the impact of the tracer–monomer interaction strengths $ε$_mer is larger than the tracer–cross-linker interaction parameter $ε$_x, which is still valid for the permeability, as obvious in Fig. 3. In the following, we elaborate on different qualitative scenarios of P(α) as a result of $ε$_x and $ε$_mer.

Regarding the weak interactions, e.g., β$ε$_x = β$ε$_mer = 0.01 [blue lines and symbols in Fig. 3(g)], monotonously decreasing curves P(α) are observed since D_in and K are decreasing functions of α. For tracer–polymer interaction strengths that are effectively attractive, i.e., $ε$_x,$ε$_mer > $ε$^c, the resulting permeability is not trivial because of the opposite trends of partitioning and diffusion. In the attractive regime, the membrane absorbs/traps more penetrants with increasing polymer density (higher binding site density), and thus, K is an increasing function of α. However, the diffusivity D_in (always below unity already because of steric obstruction) is even further decreased since the tracer mobility is reduced by the very same trapping. This anti-correlation of K and D is illustrated in Fig. S8 of the supplementary material and is the reason for the various qualitative scenarios of P(α).

At very strong interactions, e.g., β$ε$_x = 3 and β$ε$_mer = 1 [see the green symbols in Figs. 3(c), 3(f), and 3(i)], we find that the partitioning outvalues the diffusion for all α and P(α) is monotonously increasing. For moderate, yet attractive values, e.g., β$ε$_x = 1 and β$ε$_mer = 0.5 [see the orange symbols in Figs. 3(b), 3(e), and 3(h)], D_in and K are anti-correlated, and the permeability is near unity. In fact, P(α) exhibits extrema. For example, when β$ε$_x ≥ 1 and β$ε$_mer = 0.5 [panel (h)], the permeability shows a minimum and even changes from P/D₀ < 1 to P/D₀ > 1 (β$ε$_x ∈ {2, 3}). Since in this regime, K and D_in monotonously increase and decrease, respectively, we can deduce that the impact of the diffusivity dominates at low α (low ϕ_p), while for α → 0.5, the partitioning gains influence and P/D₀ is eventually above unity. This observation [P/D₀(α = 0.5) > 1] for attractive interaction parameters is further enhanced by the fact that the diffusion, as discussed in Sec. III C above, turns out to be less hindered by the trapping at very high polymer densities (α = 0.5).

Now we discuss scenarios in which the interaction parameters (⁠$ε$_x and $ε$_mer) have opposing qualities, i.e., one is attractive (larger than $ε$^c), while the other is repulsive (smaller than $ε$^c). In Sec. III B, we already recognized the non-monotonic behavior of K as a function of α for some specific interaction parameters (⁠$ε$_x and $ε$_mer). However, we also noticed that the extrema are not strongly pronounced and considered K as essentially neutral, i.e., K ≈ 1 for all α. In fact, these non-monotonicities in K(α) play a subordinate role for the occurrence of the extrema in P(α). Nevertheless, balancing $ε$_x and $ε$_mer in such a way that K ≈ 1 means that the permeability is controlled by the diffusion only, i.e., P ≈ D_in. For instance, this is the case when β$ε$_x = 0.01 and β$ε$_mer = 0.5 [compare blue symbols in Figs. 3(b), 3(e), and 3(h)]. The same will occur for $ε$_x = $ε$_mer = $ε$^c (not shown), for which K(α) ≈ 1. The extrema for P(α) are, as mentioned, caused predominantly by the anti-correlated quantities of K and D_in. Additionally, a weak maximization of P(α) is observed, e.g., for β$ε$_x = 0.01 and β$ε$_mer = 1. Since P(α) ≈ 1, we consider the effects of K and D on P as compensating in that case. This phenomenon has already been observed and discussed in previous computational studies.⁵⁴ 

Since we focused on one theory for the diffusion, precisely $DinBerg$ [Eq. (15)], and discussed two approaches for the partitioning, namely, K_μ [Eq. (10)] and K_B2 [Eq. (11)], we can compare two theories for P with the simulation results. The permeability described by $P=KB2DinBerg$ is depicted by dotted lines in Figs. 3(g)–3(i). The overall behavior is well-captured, but this theory becomes inaccurate for high interaction strengths. The theory that makes use of Eq. (10), which reads $P=KμDinBerg$⁠, is computationally more costly, yet slightly more accurate [see the solid lines in Figs. 3(g)–3(i)]. The largest deviations from the simulation results at α = 0.5 are due to the failure of $DinBerg$ at this cross-link ratio (be reminded that simulation results for α = 0.5 were not considered for the fit). The reader is referred to Secs. III B and III C for the detailed evaluation of the theories.

In this work, we investigated the effect of the cross-link ratio α on the polymer volume fraction ϕ_p and the solute permeability P of regular polymer networks by means of coarse-grained computer simulations. We particularly studied the role of the parameters $ε$_mer and $ε$_x, accounting for the specific interactions of the solutes with the chain-monomers and cross-linkers, respectively.

In the coarse-grained Langevin simulations, we imposed harmonic bond and angle potentials for bonded beads and employed the Lennard-Jones potential for non-bonded interactions in the good solvent regime (β$ε$ = 0.1). The polymer networks have been arranged as slab-like membranes in order to study the tracer permeability in accordance with the solution–diffusion model P = KD_in, i.e., by measuring the equilibrium partitioning K and the diffusivity D_in inside the polymer membrane. The permeability was studied at different cross-link ratios α as well as membrane-specific interaction parameters, namely, $ε$_mer and $ε$_x. With the chosen parameter range, the permeability P/D₀ can be tuned within two orders of magnitude approximately ranging from 10⁻¹ to 10¹.

A simple estimate of β$ε$ from all-atom simulations of our group reveals that the chosen range is appropriate for the interactions between small apolar molecules, such as butane, hexane, and benzene, with poly(N-isopropylacrylamide) (PNIPAM).⁵⁵ With typical interaction length scales of σ ≈ 0.4 nm, the proposed solute–polymer specific adsorption coefficients in that work coincide with the second virial coefficients for the LJ potential as Γ^* ≈ − 2B₂ (without distinguishing between chains and cross-linkers in this comparison). Moreover, the predictions of different qualitative scenarios for the solute partitioning in regular networks as a function of the cross-link ratio were successfully reproduced by the coarse-grained simulations herein.

In addition to the partitioning, we have presented the results for the diffusion and permeability as functions of the cross-link ratio (or like-wise as a function of the polymer volume fraction) in the present work. We found different qualitative scenarios in dependence on the interaction parameters. For given interaction strengths between the solute and the polymer, the cross-link ratios can be used to fine-tune K, D_in, and P according to the needs of a membrane device.

For a better understanding of the underlying principles, we presented scaling theories, expressed in terms of $ε$_x, $ε$_mer, α, and ϕ_p, which are in good agreement with the simulation results. These theories are based on penetrant trapping and escape processes inside the polymer network, and due to the generic approaches, they are applicable to regular extended networks with other architectures. While the two theoretical expressions provided for the partitioning stem from rather common approaches, the analytical finding for the diffusion [Eq. (15)] yields a novel, convenient expression that resolves the contributions of the component-specific interactions (⁠$ε$_mer and $ε$_x). This implies that if the relationship between α and ϕ_p for a given regular and non-collapsed network is known, a fit with the provided theory to experiments allows an interpolation and extrapolation of D_in(α) for weak and moderate interactions.

As mentioned, it is important to know how ϕ_p scales with α. In arbitrary regular networks, this is governed by the cross-linker functionality f, and we thus presented how f and α relate to $ϕpstatic$ (the polymer volume fraction of static and fully extended networks). These geometrical calculations with differing f (diamond f = 4, simple cubic f = 6, and body-centered cubic f = 8) show an increase in $ϕpstatic$ with the increase in functionality and cross-link ratio. Interestingly, we found that the relation $ϕpstatic(α)$ shows an almost universal scaling irrespective of the functionality, which is $ϕpstatic∝α2$ for small α and a trend toward linearity for higher cross-link ratios α → 0.5.

The provided theory not only yields the potential for different topologies but also can be employed to evaluate polymers with heterogeneous chains, e.g., copolymers, since the extension of the model is straightforward.

In contrast to swollen/extended networks, the physics controlling the solute dynamics in collapsed gels, amorphous polymer melts, or networks with mesh sizes comparable or smaller than the tracer size is fundamentally different. In such dense media, the penetrants are rather trapped in local cavities of the network, and the passage to another cavity is usually facilitated by a stochastic opening and closing of pathway pores.^{40,71,100–102} Such hopping diffusion strongly depends on the flexibility of the polymer and on the size and shape of the penetrant.^{40,71,102–105} Additionally, further investigation is required to evaluate how other architectures, the flexibility of the polymer, chain conformations (curving, coiling, zigzag, etc.), defects, and irregularities in the network structure compare to the presented results. Promisingly, we find substantial similarities with our results for polydisperse and more flexible networks with comparable polymer volume fractions,⁵⁴ suggesting that the semi-analytical theories of the present work are possibly transferable to polymeric networks with a less regular architecture and conformation.

Overall, this work encourages the fabrication and examination of polymer devices with regular network architectures and with an appropriate choice of the chemical composition (chain-monomers and cross-linkers) since precise control over the solute permeability and many other material properties can be achieved by tuning the cross-link ratio.

See the supplementary material for the justification of the chosen tracer concentration, details on the diffusion measurement in simulations, the calculation of the polymer volume fraction, the illustrations of possible extrema of the partitioning, the simulation results plotted against the polymer volume fraction, the simulation results with focus on small cross-link ratios, the detailed derivation of the analytical expression for the diffusion, fitting results of two theories (diffusion) at fixed cross-link ratios, and partitioning–diffusion correlation diagrams.

This project received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 646659). W.K.K. acknowledges support by a KIAS Individual Grant (No. CG076001) at the Korea Institute for Advanced Study. M.K. acknowledges financial support from the Slovenian Research Agency (Contract Nos. P1-0055 and J1-1701). This work was supported by the Deutsche Forschungsgemeinschaft via the Research Unit FOR 5099 “Reducing complexity of nonequilibrium systems.” The authors acknowledge support from the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through Grant No. INST 39/963-1 FUGG (bwForCluster NEMO).

The data that support the findings of this study are available from the corresponding author upon reasonable request.