Hairy Gels: A Computational Study

Abstract

We present results of MD and MC simulations of the equilibrium properties of swelling gels with comb-like or bottlebrush subchains and compare them to scaling-theory predictions. In accordance with theory, the simulation results demonstrate that swelling coefficient of the gel increases as a function of the polymerization degree of the main chains and exhibits a very weak maximum (or is virtually constant) as a function of the polymerization degree and grafting density of side chains. The bulk osmotic modulus passes through a shallow minimum as the polymerization degree of the side chains increases. This minimum is attributed to the onset of overlap of side chains belonging to different bottlebrush strands in the swollen gel.

1. Introduction

Brush-like macromolecules (molecular brushes) have been extensively studied theoretically and experimentally for a number of decades [1,2,3,4,5,6,7,8,9,10,11]. Both intra- and intermolecular interactions between side chains densely attached to the molecular backbone determine specific conformational and dynamic properties of molecular brushes, as compared to linear analogues. Chemical (covalent) cross-linking of brush-like polymers with subsequent swelling in a good solvent gives rise to the so called “hairy” gels with strands constituted by molecular brushes. The swelling ratios and osmotic moduli of such gels depend in a complex way on the grafting density and polymerization degree of the side chains decorating the network strands [12,13,14,15,16]. Similar structures (“hairy mesogels”) arise upon self-assembly of triblock copolymers with a comb-like or bottlebrush central block and associated terminal blocks [12,17]. In these, physically cross-linked networks the strands are formed by central blocks of copolymer (molecular brushes) that connect neighboring domains of associated blocks (network cross-links).

A molecular brush consists of a linear chain backbone with multiple side chains tethered to it. Three major synthetic approaches: (i) “grafting to” (pre-synthesized side chains are covalently attached to the backbone); (ii) “grafting through” (polymerization of the so-called macromonomers), and (iii) “grafting from” (side chains are polymerized from the backbone as macroinitiator) produced a wide variety of molecular brushes with linear grafts. An additional structural complexity can be introduced through selective gradients in grafting density or block copolymers as side chains [2,18,19].

The possibility to vary architectural parameters of molecular brushes (such as side chain length and grafting density) allows for control and adjustment of static and dynamic properties of both bulk materials and gels thereof [2,6]. Furthermore, branched architectures with bottlebrush motifs, such as barbwire [20,21,22] or dendronized polymers [23,24], have been synthesized and explored theoretically and experimentally [25,26,27]. Combining branched architecture with temperature-, pH-, and light-responsive functions in the main and side chains of molecular brushes opens new opportunities for smart material design [17,28,29,30,31,32,33,34,35]. Since properties of the constituent macromolecules play a governing role in functions of the materials, analytical theories [36,37,38,39,40,41,42] and self-consistent field numerical [43] approaches; and coarse-grained computer simulations [44,45,46,47,48,49,50,51,52,53] of branched polymers have been used to corroborate the relationships between macromolecular architecture and experimentally accessible properties of melts, solutions, and thin films comprising molecular brushes.

It was demonstrated [31] that thermo-responsive triblock copolymers with linear (L) terminal and bottlebrush (B) central blocks can produce hydrogels upon association of L-blocks in spherical domains physically cross-linking B-strands. For example, PNIPAM-bbPEG-PNIPAM triblocks self-assembled upon reaching their lower critical solution temperature (LCST) to produce two types of polymer networks: injectable hydrogels at body temperature and elastomers after water evaporation [31]. The gelation process was attributed to LCST-triggered microphase separation of the PNIPAM L-blocks, and the forming network in both hydrogels and elastomers was homogeneous, in contrast to microphase-separated linear counterparts.

Self-assembly of bottlebrush block copolymers in solutions remains a topic of intensive experimental [54,55,56,57,58,59,60,61,62,63,64] and theoretical [42,65,66] research. It was demonstrated that branched architecture of soluble blocks leads to a variety of self-assembling aggregates that could potentially serve as precursors of gels with bottlebrush strands.

The goal of this paper is to compare the recent theoretical predictions on “hairy” gels with the results of MD and MC computer simulations. In Section 2, we briefly review the scaling model of hairy gels, and summarize the theoretical predictions on its swelling behavior and mechanical properties. We then compare the theoretical predictions with the data from Monte Carlo (MC) and molecular dynamic (MD) simulations. In Section 4, we present the details of implemented MC and MD methods. In Section 3, we formulate conclusions and outline perspectives for further development in theory and modeling of bottlebrush architectures.

2. Results and Discussion

2.1. “Hairy” Polymer Gel: Scaling Model

Swelling of chemically and/or physically cross-linked networks with brush-like strands in a good solvent yields “hairy” gels. Each strand in a hairy gel constitutes a molecular brush with degree of polymerization (DP) M in its flexible backbone (main chain), and equally flexible spacers and side chains with DPs m and n, respectively (see Figure 1a). In contrast to networks with linear strands, the swelling and elastic properties of hairy gels depend not only on the cross-linking density defined by the strand DP $N=M(1+n/m)$, but also on n and grafting density ($1/m$) of side chains at a given M. If $n>m$, side chains overlap and stretch normally to the backbone due to monomer–monomer interactions, giving rise to a bottlebrush (molecular brush). If $n<m$, the side chains exhibit coil conformations, making the polymer comb-like.

A molecular brush is envisioned as a wormlike chain with thickness D (which is end-to-end distance of the side chains normally to the backbone), spacer end-to-end distance h, effective contour length $L\cong Mh/m$, and persistence length $l_p\simeq D$. In the framework of scaling model [36,42], the strand thickness D and spacer end-to-end distance h are specified by the balance between the elastic stretching of side chains and spacers, and repulsive monomer-monomer interactions in the cylindrical layer around the main chain, to give

$$D/a\cong \left\{\begin{array}{cc}{n}^{\nu }{(n/m)}^{\nu (1-\nu )/2},\hfill & m\ge m^{*}\hfill \\ n^{2\nu /(1+\nu )}{m}^{-(1-\nu )/(1+\nu )},\hfill & m\le m^{*}\hfill \end{array}\right.$$

(1)

$$h/a\cong \left\{\begin{array}{cc}{m}^{\nu }{(n/m)}^{\nu (1-\nu )/2},\hfill & m\ge m^{*}\hfill \\ m,\hfill & m\le m^{*}\hfill \end{array}\right.$$

(2)

The first and the second lines in Equations (1) and (2) correspond to partial and close to maximal stretching of spacers in the main chain, respectively (with $m^{*}\cong n^{\nu /(2+\nu )}$ separating two elasticity regimes for the backbone); a is the length of monomer unit in the backbone and side chain; and $\nu$ is Flory exponent ($\nu \approx 3/5$ in good solvent, and $\nu =1/2$ in theta-solvent). Notably, the scaling relations in Equations (1) and (2) have asymptotic character—that is, they apply only for $n\gg 1$, $n\gg m$, and $M\gg M^{*}$ with

$$M^{*}\simeq \left\{\begin{array}{cc}{m}^{1-\nu }{n}^{\nu },\hfill & m\ge m^{*}\hfill \\ n^{2\nu /(1+\nu )}{m}^{-(1-\nu )/(1+\nu )},\hfill & m\le m^{*}\hfill \end{array}\right.$$

In scaling terms, $M\simeq M^{*}$ (or, equivalently $L\simeq D$) separates bottlebrushes with $L\gg D$ from starlike polymers with $L\ll D$. In realistic experimental and simulation systems, $n\le {10}^2$ and $m\simeq 1$. Due to the local cylindrical symmetry, a strong overlap of moderately long side chains occurs only close to the backbone, and therefore the scaling dependences predicted by Equations (1) and (2) are only approached in the currently attainable range of n, $m,$ and M. However, due to still noticeable stretching of the side chains normally to the backbone ($D>a{n}^{\nu }$), a bottlebrush molecule behaves as a self-avoiding chain composed of $L/D$ impermeable subunits (“superblobs” with size D each), and its overall size $R_{eq}$ in dilute solution scales as

$$R_{eq}\cong D{\left(\frac{L}{D}\right)}^{3/5}={\left(\frac{Mh}{ma}\right)}^{3/5}{\left(\frac{D}{a}\right)}^{2/5}$$

(3)

Distribution of polymer’s density within a swollen hairy gel could be inhomogenuous. Two structural regimes, hollow mesh and filled mesh gel, are distinguished depending on the ratio between the mesh size, $R_{mesh}$, and the superblob size, D [15,16]. In the hollow mesh regime (Figure 1b), each strand constitutes a molecular brush with thickness $D\ll R_{mesh}$, and weak overlap of neighboring strands ensures hollow space (mesh) between cross-links. In this case, the mesh size $R_{mesh}$ can be evaluated using the $c^{*}$-theorem of de Gennes, [67] to give $R_{mesh}\simeq R_{eq}$. In the filled-mesh regime (Figure 1c), the neighboring strands strongly overlap, giving rise to a semi-dilute solution of side chains with an almost uniform concentration $c\simeq N{a}^3/R_{mesh}^3$. The interior of the gel is envisioned as a closely-packed array of the concentration blobs with size $\xi \left(c\right)\sim a{c}^{-\nu /(3\nu -1)}$$<D$. In this case, the equilibrium mesh size $R_{mesh}$ results from the balance between gel osmotic pressure, $\pi /k_{B}T\cong {\xi }^{-3}\left(c\right)$, and the conformational elasticity of the backbones (renormalized according to the polymer concentration c inside the gel).

To facilitate comparison between the theoretical predictions and the simulation data, we present in Figure 2 the scaling-type diagram of states for hairy gels [16] in $n,m$ log–log coordinates (with $\nu =3/5$). In the hollow-mesh regimes $C_a^{*}$ and $C_b^{*}$, each spacer in the main chain is, respectively, partially or almost fully stretched. The boundary $C_a^{*}-C_b^{*}$ corresponds thereby to $m=$$m^{*}$, indicating the onset of spacer strong stretching in individual strands at given $n.$ Similar situation occurs in the filled mesh regimes $C_a^{**}$ and $C_b^{**}$. That is, the boundary $C_a^{**}-C_b^{**}$ corresponds to the onset of spacer strong stretching in the filled-mesh regimes (semi-dilute solutions of the side chains). At the boundaries $C_b^{*}-C_b^{**}$ and $C_a^{*}-C_a^{**}$ (red lines in Figure 2), DP M of the backbone becomes equal (in scaling terms) to $M^{*}$, separating bottlebrush and starlike conformations of strands. At these boundaries (indicated by $M^{*}$), each strand comprises on the order of one superblob with size D ($L\simeq D$), as schematically shown in Figure 2.

The scaling expressions for the equilibrium mesh size $R_{mesh}/a{M}^{\nu }$, gel swelling coefficient (that is, ratio of volumes V in the swollen and dry states)

$$Q=\frac{V}{V_{dry}}\simeq \frac{R_{mesh}^3}{a^{3}N}=\frac{R_{mesh}^3}{a^{3}M}·\left(\frac{m}{n}\right),$$

(4)

with $N=M(1+n/m)$ and $n\gg m$, and osmotic bulk modulus

$$\frac{G}{k_{B}T}=c{\left(\frac{\partial \pi }{\partial c}\right)}_{c=c_{eq}}\simeq \left\{\begin{array}{cc}{R}_{mesh}^{-3}& \mathrm{hollow}\mathrm{mesh}\mathrm{regimes}{C}_a^{*},C_b^{*}\\ \xi {\left(c\right)}^{-3}& \mathrm{filled}\mathrm{mesh}\mathrm{regimes}{C}_a^{**},C_b^{**}\end{array}\right.$$

(5)

are collected in

Table 1. Asymptotic power-law dependencies for normalized mesh size $R_{mesh}$, swelling ratio Q, and bulk osmotic modulus G in a free-swelling hairy gel with DP M of strands in good solvent ($\nu =3/5$).

	${\mathit{R}}_{\mathit{mesh}}/{\mathit{aM}}^{3/5}$	${\mathit{QM}}^{-4/5}$	${\mathit{Ga}}^3{\mathit{M}}^{9/5}/{\mathit{k}}_{\mathit{B}}\mathit{T}$
$C_a^{*}$	${(n/m)}^{9/25}$	${(n/m)}^{2/25}$	${(n/m)}^{-27/25}$
$C_a^{**}$	${(n/m)}^{9/25}$	${(n/m)}^{2/25}$	${(n/m)}^{-9/50}$
$C_b^{*}$	${(n/m)}^{3/10}{m}^{1/5}$	${(n/m)}^{-1/10}{m}^{3/5}$	${(n/m)}^{-9/10}{m}^{-3/5}$
$C_b^{**}$	$M^{2/5}$	$M^{6/5}{(n/m)}^{-1}$	$M^{-27/10}{(n/m)}^{9/4}$

with corresponding exponents specified for $\nu =3/5$.

As follows from Table 1, an increase in DP n of the side chains at fixed DP M of the backbone leads to the monotonous increase in mesh size, $R_{mesh}\sim n^{\beta }$, in regimes $C_a^{*}$ and $C_a^{**}$ (exponent $\beta =9/25=0.36)$, $C_b^{*}$ (exponent $\beta =3/10$), approaching full extension, $R_{mesh}\sim aM$, in regime $C_b^{**}$ (exponent $\beta =0$). At the same time, both swelling ratio Q and osmotic bulk modulus G exhibit non-monotonic dependences on n. Swelling coefficient $Q\left(n\right)\sim n^{3\beta -1}$ passes through a maximum at the boundaries $C_a^{*}-C_b^{*}$ and $C_a^{**}-C_b^{**}$, and osmotic modulus $G\left(n\right)$ passes through a minimum upon crossing the boundaries of regime $C_b^{**}$. The maximum in $Q\left(n\right)$ dependence is weak due to a small value of the exponent which changes from $+0.08$ to $-0.1$ at the $C_a^{*}-C_b^{*}$ boundary, and this maximum could even disappear if the apparent exponent ${\beta }_{app}<1/3$. Recall that the values of exponents in Table 1 were calculated in the limit of $n/m\gg 1$. The predicted sharp (jumpwise) minimum in $G\left(n\right)$ at the boundaries of regime $C_b^{**}$ could be smoothed in computer simulations, leading to shift in the minimum location far left (to smaller values of n).

2.2. Computer Simulations of Hairy Gel

We used both MC and MD simulations to explore the equilibrium-swelling behavior of hairy gels in a good solvent ($\nu =3/5$). In Figure 3, we present the simulation box which has the same geometry in MC and MD simulations. That is, 16 polymer chains were connected to a diamond-like network by eight tetrafunctional cross-linking units, and the network element was put into cubic simulation box of the volume V with periodic boundary conditions to emulate an infinite polymer network. It is worth mentioning that MD and MC simulations were performed in $NVT$ and $NPT$ (with pressure $P=0$ and V fluctuating) ensembles, respectively.

In Figure 4a,b, we present the scaling-type diagrams with positions of the boundaries calculated specifically for DP $M=37$ and DP $M=13$ of the backbone, respectively, together with set of symbols $\{n,m\}$ marking the architectural parameters of bottlebrush strands modeled in MC simulations. In Figure 4c, a similar diagram is presented for $M=30$ with set of symbols $\{n,m\}$ marking the parameters of strands modeled in MD simulations.

While positions of the boundaries in Figure 4 are specified with accuracy of the numerical prefactors on the order of unity, it is still expected that MC simulations for $M=37$ (Figure 4a) with small values of $1\le m\le 4$ and increasing n (up to $n=128$) cover regimes $C_a^{*}$ and $C_b^{*}$. In contrast, Figure 4b indicates that for $M=13$, MC data pass through regimes $C_a^{*}$$C_b^{*}$, and even enters regime $C_b^{**}$. According to Figure 4c, MD simulations cover a considerably larger interval of m-values (up to $m=50$). However, many of the symbols correspond to comb-like strands (area with $n/m<1,$ shaded gray), and the rest of the data probes only regime $C_a^{*}$.

2.2.1. Average End-to-End Distances of Strands and Side Chains

In Figure 5, we plot the average mesh size $R_{mesh}\simeq {\langle V\rangle }^{1/3}$ of a free-swelling gel, and the average end-to-end distance $\langle r\rangle \simeq r_{sc}\simeq D$ of the side chains obtained in MC simulations as a function of n for $m=1,2,3,$ and 4 in log–log coordinates.

At $n\gtrsim 10$ (when the side chains stretch normally to the backbone), different m-values produce almost parallel lines in $R_{mesh}\left(n\right)\sim n^{\beta }$ dependences, with apparent slopes ${\beta }_{app}=0.22-0.25$, i.e., smaller than the theoretically predicted $\beta$ (i.e., $0.36$ in regime $C_a^{*}$ and $0.30$ in regime $C_b^{*}$). For $r_{sc}\left(n\right)\sim n^{\alpha }$, the correspondence between the theoretical exponents $\alpha$ (i.e., $0.75$ in regimes $C_b^{*}$ and $0.72$ in regime $C_a^{*}$) and apparent slopes is expectably [50] worse, consistent with smaller apparent values of ${\beta }_{app}$. Notably, at $M=37$, $R_{mesh}\left(n\right)$, and $r_{sc}\left(n\right)$ do not intersect at fixed m at any considered n (see Figure 5a), indicating that the hairy gel remains in the hollow mesh state. In contrast, at $M=13$, $R_{mesh}\left(n\right)$ and $r_{sc}\left(n\right)$ become close at the largest values of n (see Figure 5b), indicating that the hairy gel approaches the filled-mesh state. In Figure 5c, we plot normalized mesh size $R_{mesh}/a{M}^{3/5}$ as a function of ($n/m)$ for $M=13,19,25$, and 37 to demonstrate how MC data collapse on mastercurve with slope ${\beta }_{app}=0.25$, with maximal deviations for the smallest value of $M=13$ (presumably approaching regime $C_b^{**}$ with $R_{mesh}/a{M}^{3/5}\sim {(n/m)}^0$).

In Figure 6, the equilibrium strand end-to-end distance $R_{mesh}\left(n\right)\sim <V{>}^{1/3}$ is presented for a series of m-values and fixed $M=30,$ as obtained from MD simulations. While values of $n<m$ correspond to comb-like strands, the data for relatively small m and $n>10$ are used in Figure 6 to evaluate apparent exponent ${\beta }_{app}$ in the dependence $R_{mesh}\left(n\right)M^{-3/5}\sim {(n/m)}^{{\beta }_{app}}$ to give ${\beta }_{app}\approx 0.26$, in accordance with the results of MC simulations.

2.2.2. Gel Swelling Coefficient

In Figure 7, we present the normalized equilibrium swelling coefficient $Q{M}^{-4/5}=V{M}^{-4/5}/a^{3}N$ as a function of $n/m$ predicted by the scaling model (Figure 7a) and obtained from MC simulations (Figure 7b) in log–log coordinates. The power-law dependences $Q{M}^{-4/5}$ in Figure 7a were calculated for $M=37$ with slopes indicated in Table 1 (that is, $2/25,-1/10,$ and $-1$ in regimes $C_a^{*},C_b^{*},C_b^{**}$, respectively) and $m=1,2,3$, with all numerical prefactors assigned unity. For $m=2$, Q is shown by red lines; for $m\ge 4$ only regimes $C_a^{*},C_b^{**}$ are feasible. The maximum predicted for $M=37$ corresponds to $Q{M}^{-4/5}$$\approx 1.38$. A decrease in M shifts location of the maximum to the left and makes it less pronounced.

MC simulations data in Figure 7b show the dependence of $Q{M}^{-4/5}$ on $n/m$ for series of M-values. The predicted weak maximum is not well pronounced; mostly, its decreasing (right) branch is seen in MC simulations. As strands with $M=37$ are not expected to enter regime $C_b^{**}$ (see Figure 4a), the slope $dQ/dn$ is far from the predicted exponent $-1$. However, for smaller $M=13$ for which regime $C_b^{**}$ is feasible (see Figure 4b), the slope of MC data approaches $-1$, as predicted.

The data from MD simulations in Figure 8 also indicate weak dependence of swelling coefficient Q on $n/m$. Here, the normalized swelling coefficient $Q{M}^{-4/5}$ is presented as a function of $n/m$ for a wider variety of m-values and two values of $M=30$ (red symbols) and $M=100$ (violet symbols), and probing regime $C_a^{*}$ in which the predicted slope is $2/25=0.08.$ For $M=30$, the data in Figure 8 with different m-values remain rather scattered; the increase in backbone DP M up to $M=100$ decreases the scattering of the data that collapse on mastercurve with close to zero slope. Notably, the numerical prefactor in $Q{M}^{-4/5}$ versus $n/m$ dependences is close to unity in both MC and MD simulations, consistent with the scaling model.

2.2.3. Osmotic Bulk Modulus

In Figure 9a we present the MC data for the gel osmotic bulk modulus, $G=c{(\partial \pi /\partial c)}_{c=c_{eq}}$. Here, $\pi \left(c\right)$ is osmotic pressure, and $c_{eq}\simeq Q^{-1}$ is the equilibrium concentration of monomer units in the hairy gel. As is seen in Figure 9a, the osmotic modulus G decreases with increasing DP M of the strand backbone, in qualitative agreement with the theoretical predictions in Table 1. However, the predicted minimum in $G\left(n\right)$ dependence is not detected in MC simulations. An increase in M flattens $G\left(n\right)$ dependence, only pointing at the possibility of minimum formation.

A conjecture is that the sharp minimum predicted by the scaling model could be smoothed in MC simulations, as schematically illustrated by the dashed lines in Figure 9b. Another factor is the effect of perturbed dense cylindrical layers circumventing the backbones. Although Figure 9b is merely schematic and does not specify the shape or position of the dotted lines, it is clear that the theoretically predicted exponents for $G\left(n\right)$ would not work in the smoothed region, and therefore, collapse of G-data on the mastercurve is not expected. However, this schematic indicates that minimum in $G\left(n\right)$ dependence can move left and fall out of the considered ranges of the strand architectural parameters.

3. Conclusions

Based on the scaling description of molecular brushes, the theory specifies power-law dependences for the equilibrium properties of “hairy” gels with bottlebrush strands. Four structural regimes of hairy gels were distinguished: two hollow-mesh regimes, $C_a^{*}$ and $C_b^{*}$, with partially and almost fully stretched spacers; and two filled mesh regimes, $C_a^{**}$ and $C_b^{**}$, with partially or almost fully stretched backbones embedded in a semi-dilute solution of side chains. Two supplementary computer simulation techniques, MC and MD, were used to probe the structure of free-swelling in good solvent gel with varied strand architectural parameters $\{M,n,m\}$ to corroborate the theoretical model. It was demonstrated that MC simulations of hairy gels with short spacers ($m\le 4$) could cover regimes $C_a^{*}$ and $C_b^{*}$, and approach regime $C_b^{**}$. MD simulations of gels with a wider variety of m-values could probe regime $C_a^{*}$. In all cases, limited extension of the gel regimes (Figure 4) and restriction $n/m\gg 1$ make it challenging to compare the predicted asymptotic values of exponents (Table 1) with apparent exponents from MC and MD simulations. In addition, it is not surprising that apparent exponents deviate from the theoretical ones for the considered set of strand parameters $\{M,n,m\}$.

According to the scaling model, the gel mesh size $R_{mesh}\sim M^{\nu }{n}^{\beta }$ increases monotonously with increasing of both the DP n of the side chain and the DP M of the strand backbone. The apparent exponent ${\beta }_{app}\approx 0.25-0.26$ estimated from MC and MD simulations was smaller but reasonably close to the theoretical values, $\beta =0.30-0.36$. Exponent $\nu$ was close to $3/5$ as expected for good solvent conditions. The thickness D of the bottlebrush strand increased with n. However, an approach to the theoretical exponent $3/4$ requires $n\gtrsim {10}^3$ (i.e., much larger n-values than have been implemented in MC simulations), and the correspondence between apparent and scaling exponents in $D\left(n\right)$ dependence was poor.

The theory predicts that the swelling ratio Q of the hairy gels is controlled primarily by the DP M of the strand backbone, whereas the dependence on the DP n and grafting density $1/m$ of the side chains is weak. This prediction is in agreement with the MC and MD simulation data. The theory also predicts non-monotonic dependences of osmotic modulus G and swelling ratio Q of hairy gels on DP n of the side chains. The osmotic modulus G passes through a minimum corresponding to the overlap threshold of the side chains emanating from different strands. In contrast, swelling ratio Q passes through a maximum upon onset of strong stretching of spacers in nonoverlaping strands. The MC computer simulations did not find pronounced extrema in $Q(n/m)$ and $G(n/m)$ dependences and demonstrated only the decreasing branch in Q-dependence (with slope approaching $-1$, predicted in regime $C_b^{**}$) and the increasing branch in G-dependence (with slope approaching $9/4$, predicted in regime $C_b^{**}$). Lack of maximum in $Q\left(n\right)$ dependence can be explained by the values of apparent exponent ${\beta }_{app}<1/3$ in $R_{mesh}\sim n^{{\beta }_{app}}$ dependence that automatically eliminates the increasing branch in $Q=R_{mesh}^3/a^{3}N$$\sim n^{3{\beta }_{app}-1}$. The relatively small values of ${\beta }_{app}$ could follow from the limited extension of the scaling regimes for currently considered set of $\{M,n,m\}$ parameters. A wider set of architectural parameters $\{M,n,m\}$ is desirable to confront the scaling model of hairy gel with better accuracy.

4. Materials and Methods

The model of the gel as a network of 16 linear polymer chains, each consisting of M monomer units. These polymer chains are connected to a diamond-like network by eight cross-linking units. The side chains of the lenght n are grafted on each m-th monomer of the main chain (see Figure 3). The network was put into cubic simulation box of the volume $V_{gel}$ with periodic boundary conditions, which virtually emulates an infinite polymer network.

4.1. Monte Carlo Simulations

In Monte Carlo simulations, we used a variant of the Hamiltonian (originally called hybrid [68]) Monte Carlo (HMC) method and coarse-grained (CG) models. The HMC uses Hamiltonian dynamics to sample probability distribution $exp(-H/k_{B}T)$, where $H(p,q)$ is the Hamiltonian, q are generalized coordinates, and p are generalized momenta. In the simplest case, a separable Hamiltonian $H(p,q)=K\left(q\right)+V\left(p\right)$ is used, where K is kinetic energy and V is potential energy for sampling Boltzmann distribution $exp(-V\left(q\right)/k_{B}T)$. The momenta p can be sampled directly, and Hamiltonian dynamics $\dot{q}=\partial H/\partial p$, $\dot{p}=-\partial H/\partial q$ are followed for some time to prepare a new proposal for a Metropolis step [69] with acceptance probability $P_{\mathrm{acc}}=min(1,exp(-\Delta H/k_{B}T))$. For the exact dynamics $\Delta H$, it would be zero (time independent Hamiltonian is conserved), and all proposals would be accepted. In practice, instead of exact dynamics, we must use an approximate numerical evolution with numerical integrators. Fortunately, time-reversible, phase-space volume-preserving integrators are available, and their procedure is exact, and any bias due to approximate dynamics is then removed by Metropolis rejections. Moreover, the Hamiltonian for dynamic evolution can be different from the targeted one. This allows for great variability of the method. In the case of polymer simulations, a big improvement can be achieved by a suitable transformation of variables. The potential energy for the standard Rouse (harmonic) model is diagonal in bond vector coordinates $r_i-r_{i-i}$. With the freedom of choosing the evolution Hamiltonian, we use new canonically conjugated momenta with bond vector coordinates, and it allows for sampling all normal modes at the same rate [70]. For other potentials and approximate integrators, the method somewhat deteriorates, but it is still much better than simple molecular dynamics and leads to a smaller scaling exponent of the integrated autocorrelation time of the end-to-end distance with chain length, and thus is arbitrarily faster.

In our CG model we use two types of interactions, bonding and non-bonding. For the bonding one, we use a variant of finitely extensible nonlinear elastic (FENE) potential [71] in a form $u_{\mathrm{F}}\left(r\right)=-e{k}_{\mathrm{B}}Tlog((r_0+d-x)(x-r_0+d)/d^2)$ for $r\in (r_0-d,r_0+d)$ and infinite elsewhere, with mild singularities at $r_0\pm d$. It is a relatively small complication for simple molecular dynamics where small time-steps are used that potential does not exist outside $(r_0-d,r_0+d)$. In HMC, we use the advantage of much longer time-steps that could lead to stepping out of the definition interval. We resolve this in the spirit of HMC and prepare a new potential finite, everywhere being $u_{\mathrm{F}}\left(r\right)$ on $(r_0-d+ϵ,r_0+d-ϵ)$ and $v\left(r\right)+u(r_0+d-ϵ)-v(r_0+d-ϵ)$, where $v\left(r\right)=e{k}_{\mathrm{B}}T{(r_0/d)}^2{(r/r_0-1)}^2$ elsewhere that is used for dynamic evolution in HMC and the original $u_{\mathrm{F}}\left(r\right)$ is used for the Metropolis step. By this combination, we achieve the exact sampling with $u_{\mathrm{F}}\left(r\right)$ while avoiding numerical problems. For non-bonding potential, we use soft repulsive potential $u_S=e{((r-c)/r)}^2$ for $r\in (0,c)$ and 0 elsewhere, where $e>0$ is an energy parameter and c is a cutoff. The potential is smooth at cutoff. Frequently, the standard Lennard–Jones potential [72], originally suggested for simulations of neon, is used for this purpose. Although the $1/r^6$ part can be justified for two atoms (with at least one in an S-state) at very long distances, the $1/r^{12}$ part lacks such a simple justification, the less for the potential that should represent CG macromolecular system in a solution. By taking a systematically coarse-grained model based on all-atom empirical force-field and fitting the corresponding part for small distances, we find a much smaller exponent [73] of about 2, as for our $u\left(r\right)$.

The simulations were performed in $NPT$ ensemble by adding volume-changing moves [74]. The simulations started from extended conformations of gel backbone and coiled side-chains as corresponding to a good solvent. A thorough equlibration was performed for box sizes, energies, and end-to-end distances. The free (tuning) parameters of our HMC method were roughly set up according to optimal acceptance probability and the smallest integrated autocorrelation times. After equilibration (as illustrated by Figure 10), samples were accumulated, and standard deviations of their averages were estimated using the blocking method [75].

4.2. Molecular Dynamics Simulations

Each pair of the particles interact via the truncated Lennard–Jones interaction potential, which imposes strong repulsion between all particles at short distances:

$$\begin{array}{c}\hfill U_{\mathrm{LJ}}\left(r\right)=\left\{\begin{array}{cc}4\epsilon \left({\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right)\hfill & \mathrm{if}r<r_{\mathrm{cut}}\hfill \\ 0\hfill & \mathrm{elsewhere}\hfill \end{array}\right.,\end{array}$$

(6)

where r is the interparticle distance, $\sigma =0.35$ nm is a chosen characteristic size of the particles, $\epsilon =k_{B}T$ is the depth of the potential, and $r_{\mathrm{cut}}$ is the cut-off distance beyond which the potential is set zero.

The bonds connecting the gel to a network are modeled using finite extension nonlinear elastic potential (FENE):

$$U_{\mathrm{FENE}}\left(r\right)=-\frac{1}{2}K\Delta r_{\max }^{2}ln\left[1-{\left(\frac{r-r_0}{\Delta r_{\max }}\right)}^2\right],$$

(7)

where r is the distance between the bonded segments, K is the magnitude of their interaction, $\Delta r_{\max }$ is the maximal stretching length of the bond, and $r_0$ is the equilibrium bond length. In our simulations we set $K=10kT/{\sigma }^2$, $\Delta r_{\max }=3\sigma$ and $r_0=0.0$ [76].

The Langevin thermostat [71] was used to guarantee the constant temperature of the system, $T=300$ K. The two additional terms to force in equation of motion were added.

$${\mathbf{f}}_i=-\gamma {\mathbf{v}}_i\left(t\right)+\sqrt{2\gamma kT}{\mathbf{\eta }}_i\left(t\right),$$

(8)

where the first term corresponds to constant friction, with $\gamma$ being a friction coefficient, and the second one corresponds to random thermal force, with ${\mathbf{\eta }}_i$ being a normally distributed random vector; ${\mathbf{v}}_i$—velocity of i-th particle, t—time.

In order to calculate the swelling ratio of hydrogel, we performed a series of simulations of the gel in a box of different volumes V. Each simulation results in a particular value of pressure P.

Our target observable is the free swelling equilibrium state, i.e., at the state where the applied to the gel pressure equals to zero. In order to localize this state, we first plotted $P\left(V\right)$ dependence; then, using the least-squares method, we drew a line passing through the nearest to the V axis points. The place where this line crosses the V-axis defines the volume of free swelling equilibrium state.