Thermoreversible Associating Polymer Networks: I. Interplay of Thermodynamics, Chemical Kinetics, and Polymer Physics

TL;DR

This paper introduces a simulation model for thermoreversible associating polymer networks that effectively separates thermodynamic and kinetic effects, providing insights into their dynamics and mechanical behavior.

Contribution

It presents and validates a novel simulation model capable of analyzing thermodynamics, kinetics, and polymer physics interactions in associating polymer systems.

Findings

Crossover from diffusion-limited to kinetically-limited bond recombination influences dynamics.

Model accurately captures nonequilibrium and nonlinear mechanical properties.

Polymer connectivity affects bond dynamics and material behavior.

Abstract

Hybrid molecular dynamics/Monte Carlo simulations used to study melts of unentangled, thermoreversibly associating supramolecular polymers. In this first of a series of papers, we describe and validate a model that is effective in separating the effects of thermodynamics and chemical kinetics on the dynamics and mechanics of these systems, and is extensible to arbitrarily nonequilibrium situations and nonlinear mechanical properties. We examine the model's quiescent (and heterogeneous) dynamics, nonequilibrium chemical dynamics, and mechanical properties. Many of our results may be understood in terms of the crossover from diffusion-limited to kinetically-limited sticky bond recombination, which both influences and is influenced by polymer physics, i. e. the connectivity of the parent chains.

Figures

Tables

Table 1. Table 1: Variation of τ s ​ b subscript 𝜏 𝑠 𝑏 \tau_{sb} , τ ∗ superscript 𝜏 \tau^{*} and D 𝐷 D with N 𝑁 N , T 𝑇 T , h ℎ h and τ M ​ C subscript 𝜏 𝑀 𝐶 \tau_{MC} . Times are in units of τ L ​ J subscript 𝜏 𝐿 𝐽 \tau_{LJ} . Statistical errors are roughly ± 2 % plus-or-minus percent 2 \pm 2\% or less. All systems have N c ​ h = 1400 subscript 𝑁 𝑐 ℎ 1400 N_{ch}=1400 and data are averaged over multiple statistically independent states. * denotes τ 0 = .5 ​ τ L ​ J subscript 𝜏 0 .5 subscript 𝜏 𝐿 𝐽 \tau_{0}=.5\tau_{LJ} results. − ⁣ − -- indicates calculation is prohibitive or we have insufficient data. Results for D 𝐷 D in N = 1 𝑁 1 N=1 systems are not presented because they are negligibly affected by sticky bonding.

System	$N$	$k_{B}T/u_0$	$h/k_{B}T$	${\tau }_{MC}$	${\tau }_{sb}$	${\tau }^{\ast }$	${\tau }^{\ast }/{\tau }_{sb}$	${10}^5{\tau }_{LJ}D/a^2$
A	50	1.0	10	1	$1.55\cdot {10}^3$	$7.05\cdot {10}^3$	4.5	18.3
B	50	1.0	10	10	$1.54\cdot {10}^4$	$2.63\cdot {10}^4$	1.7	8.44
C	50	1.0	10	100	$1.53\cdot {10}^5$	$1.78\cdot {10}^5$	1.2	2.17
E*	50	1.0	11.25	0.5	$2.40\cdot {10}^3$	$2.28\cdot {10}^4$	9.5	9.31
DE*	1	1.0	11.25	${0.5}^{\ast }$	$3.35\cdot {10}^3$	$2.23\cdot {10}^4$	6.7	N/A
F	50	1.0	11.25	1	$4.82\cdot {10}^3$	$2.95\cdot {10}^4$	6.1	6.89
DF	1	1.0	11.25	1	$6.69\cdot {10}^3$	$2.67\cdot {10}^4$	4.0	N/A
G	50	1.0	11.25	2	$9.78\cdot {10}^3$	$4.05\cdot {10}^4$	4.1	6.07
DG*	1	1.0	11.25	2	$1.34\cdot {10}^4$	$3.46\cdot {10}^4$	2.6	N/A
H	50	1.0	11.25	5	$2.40\cdot {10}^4$	$6.42\cdot {10}^4$	2.7	4.81
DH*	1	1.0	11.25	5	$3.36\cdot {10}^4$	$5.50\cdot {10}^4$	1.6	N/A
I	50	1.0	11.25	10	$4.73\cdot {10}^4$	$9.70\cdot {10}^4$	2.1	3.11
DI*	1	1.0	11.25	10	$6.42\cdot {10}^4$	$8.52\cdot {10}^4$	1.3	N/A
J	50	1.0	11.25	${10}^{1.5}$	$1.55\cdot {10}^5$	$2.32\cdot {10}^5$	1.5	1.95
DJ*	1	1.0	11.25	${10}^{1.5}$	$2.05\cdot {10}^5$	$2.30\cdot {10}^5$	1.1	N/A
K	50	1.0	11.25	100	$4.82\cdot {10}^5$	$6.28\cdot {10}^5$	1.3	0.94
DK*	1	1.0	11.25	100	$6.70\cdot {10}^5$	$6.85\cdot {10}^5$	$\sim 1.02$	N/A
L	50	0.6	11.25	1	$1.88\cdot {10}^3$	$8.02\cdot {10}^3$	4.3	$--$
M	50	0.6	11.25	10	$1.91\cdot {10}^4$	$3.15\cdot {10}^4$	1.6	$--$
N	50	0.6	11.25	100	$1.93\cdot {10}^5$	$2.22\cdot {10}^5$	1.2	$--$
O	50	1.0	12.5	1	$1.53\cdot {10}^4$	$1.23\cdot {10}^5$	8.0	2.79
P	50	1.0	12.5	10	$1.55\cdot {10}^5$	$3.58\cdot {10}^5$	2.3	$--$
Q	50	1.0	12.5	100	$1.53\cdot {10}^6$	$2.07\cdot {10}^6$	1.3	$--$
R	50	1.0	13.75	1	$4.90\cdot {10}^4$	$5.42\cdot {10}^5$	11	$--$
S	50	1.0	13.75	100	$4.62\cdot {10}^6$	$7.34\cdot {10}^6$	1.6	$--$
T	50	1.0	15	1	$1.56\cdot {10}^5$	$2.49\cdot {10}^6$	16	$--$
U	50	1.0	15	100	$1.48\cdot {10}^7$	$2.53\cdot {10}^7$	1.7	$--$
V	50	1.0	16.25	1	$5.0\cdot {10}^5$	$1.1\cdot {10}^7$	22	$--$

Equations

$$U_{sb}(r,h)=U_{FENE}(r)-U_{FENE}(r_0)-h,$$

$$\begin{array}[]{ccc}V-W+X&\rightarrow&V-X+W\\ V-W+X-Y&\rightarrow&V-Y+X-W\end{array}$$

$$\begin{array}[]{rcl}\left[A\right]&\equiv&\rho c_{st}(1-p_{active}),\\ &&\\ \left[A_{2}\right]&\equiv&\rho c_{st}p_{active}/2,\end{array}$$

$$\begin{array}[]{rcccl}K_{eq}&\equiv&\displaystyle\frac{[A_{2}]}{[A]^{2}}&\equiv&\displaystyle\frac{p^{*}}{2\rho c_{st}(1-p^{*})^{2}}\end{array}$$

$$\begin{array}[]{rcl}A+A&\stackrel{{\scriptstyle\overset{k_{f}}{\longrightarrow}}}{{\underset{k_{b}}{\longleftarrow}}}&A_{2}\end{array}$$

$$\begin{array}[]{rcl}k_{f}[A]^{2}&=&k_{b}[A_{2}].\end{array}$$

$$\begin{array}[]{rcl}k_{f}&=&\alpha\exp(-\delta(h)/k_{B}T),\\ &&\\ k_{b}&=&\beta\exp(-(h+\delta(h))/k_{B}T),\\ \end{array}$$

$$K_{eq}^{TS}\equiv \frac{k_f}{k_b}=\frac{\alpha }{\beta }exp(h/k_{B}T).$$

$$<(\delta \vec{r}{)}^2(t)>\equiv a_{cage}^{2}f(t).$$

$$\frac{{\tau }^{\ast }}{{\tau }_{sb}}=C+\frac{K}{({\tau }_{MC}{)}^x},$$

$${\tau }^{\ast }=\frac{C+K{\tau }_{MC}^{-x}}{k_b}.$$

$${\tau }_{deloc}={\tau }_{sb}\frac{y\log y}{y-1},$$

$$\frac{{\tau }_{deloc}}{{\tau }^{\ast }}=\frac{\log y}{y-1},$$

$$\Lambda (t)=\frac{3<\delta r^4(t)>}{5<\delta r^2(t){>}^2}-1,$$

$$\begin{array}[]{rcl}\dot{[A_{2}]}&=&k_{f}[A]^{2}-k_{b}[A_{2}],\end{array}$$

$${\dot{p}}_{active}=2k_f\rho c_{st}(1-p_{active}{)}^2-k_b{p}_{active}.$$

$$p_{active}(t)=d-\frac{\left(d^2-1\right)\tanh \left(2z\sqrt{d^2-1}t\right)+d\sqrt{d^2-1}}{d\tanh \left(2z\sqrt{d^2-1}t\right)+\sqrt{d^2-1}},$$

$$\begin{array}[]{l}\sigma_{z}(\lambda)=G_{e}exp(-ln(\lambda)/\dot{\epsilon}\tau_{sb})\times\\ \\ \left[g(\lambda)+\displaystyle\frac{1-\lambda^{2}}{1-2\dot{\epsilon}\tau_{sb}}+\displaystyle\frac{1/\lambda-1}{1+\dot{\epsilon}\tau_{sb}}\right]\end{array}$$

Full text

I Introduction

Flexible synthetic polymers have long been of fundamental scientific interest because many of their properties arise from a few universal features like the topological connectivity, random-walk like structure, and excluded volume of the chain molecules. Less universal are the various attractive, “associative” intermolecular interactions Rubinstein and Dobrynin (1999); Brunsveld et al. (2001) ranging from weak dispersion forces to strong covalent chemical bonds (in chemically crosslinked systems). Examples include hydrogen bonding, electrostatic attractions, and effective attractions driven by incompatibility with a solvent. These interactions lead to formation of supramolecular structures ranging from micelles to network gels.

Associating polymers (APs) differ from simple homopolymers in that chains contain a (typically fairly low) fraction of “sticky” monomers, which are different from the majority-species monomers. The sticky monomers form “sticky bonds” with each other via associative interactions weaker than permanent covalent bonds. These lead to formation of supramolecular aggregrates. Unlike the closely related “living” or “equilibrium” polymers, the degree of polymerization of AP “parent” chains is fixed (in time) by permanent covalent backbone bonds.

The lifetime of the sticky bonds is finite. Depending on the nature of the associative interactions and ambient conditions (e. g. temperature, concentration), the supramolecular topology may be practically permanent, in which case the system forms a “chemical gel” (i. e. a crosslinked rubber), or so short-lived that the system is indistinguishable from a simple polymer solution, melt or glass. Between these limits, when the topology of the associated supramolecular aggregates changes on a time scale comparable to the experiment, these systems form complex fluids with fascinating dynamical and mechanical properties Rotello and S. Thayumanavan (2008).

At fixed ambient conditions, the time scales for topological changes in associating polymer systems are in principle set by three independent factors: the (a) thermodynamics (i. e. energetic strength relative to $k_{B}T$), (b) the “chemical kinetics” of the sticky bonds, and (c) the underlying non-associative polymer physics. Thermodynamics set “static” quantities such as the size of the supramolecular aggregates and hence the position of the system relative to the percolative gelation transition. Kinetics set relaxation times through their effect on the rates of formation and breaking of sticky bonds. Polymer physics alters the dynamics through such effects as the random-walk-like structure and uncrossability of chains, which give rise to the systems’ underlying Rouse or reptation dynamics Doi and Edwards (1986). The interplay of (a)-(c) allows for the design of materials with exquisitely tunable rheological response. For this reason, APs have been the focus of intense experimental and theoretical study over the past two decades; see Refs. Brunsveld et al. (2001); Rotello and S. Thayumanavan (2008); Annable et al. (1993); Binder and Zirbs (2007); Tanaka (2002) for reviews.

Changes in ambient conditions lead to “thermoreversible” property changes unique to AP systems, e. g. extremely sharply decreasing viscosity upon decreasing concentration or increasing temperature. These changes can be tuned (engineered), so APs have great potential as “smart” materials Brunsveld et al. (2001); Binder and Zirbs (2007); Rotello and S. Thayumanavan (2008); de Greef and Meijer (2008); Cordier et al. (2008) in which the change of lifetime or concentration of the sticky bonds with ambient conditions leads to useful products. Applications include temperature-sensitive adhesives, coatings for heat-sensitive materials, and generally enhanced melt processability relative to conventional polymers de Greef and Meijer (2008).

Static thermodynamic properties of AP systems, in particular the percolative gelation transition and local structural changes arising from associative interactions, have been extensively studied. Analytic theories provide a good understanding of homogenous systems, and emerging numerical techniques such as self consistent field theoretic simulations Lee et al. (2007) and reaction-ensemble DPD Lisal et al. (2009) show promise for investigating inhomogeneous systems.

However, the dynamical, mechanical, and nonequilibrium properties of associating polymer gels and networks remain poorly understood. Time dependent properties obviously depend on kinetics, and in addition to being of fundamental scientific interest, a better understanding of them may prove important in developing new applications of associating polymer systems such as self healing materials Cordier et al. (2008). The situation becomes particularly complex when the lifetime of the sticky bonds is not long or short compared to the “polymeric” relaxation times; this regime has been studied rather extensively for linear equilibrium polymers (see e. g. Refs. Cates (1987); O’Shaughnessy and Yu (1995); com (a); Nyrkova and Semenov (2007); Stukalin et al. (2008)), but much less so for networks.

Analytic and quasi-analytic approaches to AP dynamics and mechanics, e. g. foo (a) Refs. Green and Tobolsky (1946); Cates (1988); Baxandall (1989); Leibler et al. (1991); Tanaka and Edwards (1992); Groot et al. (1996); Rubinstein and Semenov (1998); Vaccaro and Marruci (2000); Rubinstein and Semenov (2001); Tanaka (2002); Semenov (2006); ind ; Semenov et al. (2007); Tripathi et al. (2006), have made many useful, experimentally verifiable predictions, including nonlinear behaviors such as shear thickening and strain hardening Pellens et al. (2004a); Tripathi et al. (2006). In the general case, however, the complex interplay of sticky bond thermodynamics and kinetics with the underlying polymer physics in these systems is almost certainly beyond the reach of analytic theory. For the sake of tractability, theories have generally neglected one or more features of AP systems that are likely essential to capturing their behavior under certain ambient conditions. For example, as temperature drops towards the glass transition, attractive, non-associative interactions, such as van der Waals forces between non-sticky monomers, become increasingly important Baschnagel and Varnik (2005). Moreover, virtually all analytic treatments have thus far been restricted foo (b) to homogeneous AP systems; the majority focus on the “telechelic” case of APs with only 2 sticky monomers per parent chain (one on each end). We believe that inhomogeneous AP systems are the potentially the most interesting and useful, e. g. because inhomogeneities serve to localize sticky monomer concentration and network connectivity, which in turn can broaden the relaxation spectrum Kumar et al. (2007); Feldman et al. (2009).

The above set of potentially essential features of APs is not treated microscopically by existing theories, but can be readily captured by particle-based simulations. This is the first of a series of simulation studies, the goal of which is to elucidate the separate effects of sticky bond thermodynamics, kinetics, and other underlying polymer physics on the dynamical, mechanical, and nonequilibrium properties of associating polymers.

Previous particle based simulations of AP networks focusing on dynamical properties, e. g. foo (a) Refs.  Groot and Agterof (1994); Nguyen-Misra and Mattice (1995); Y.Liu and Pandey (1996); Milchev et al. (1999); Kumar and Douglas (2001); Ayyagari et al. (2004); Guo and Luijten (2005); Manassero et al. (2005); Loverde et al. (2005); com (b); Baljon et al. (2007); Goswami et al. (2008) have produced a wealth of interesting ideas and results. However, for various reasons outlined in Section II.2, the methods employed in these previous studies are not suitable for the full range of problems we wish to consider. A more versatile model for APs requires the combination of: i) realistic dynamics, ii) applicability to far-from-equilibrium conditions, iii) controllable sticky binding topology, iv) variable chemical kinetics and v) the ability to treat inhomogeneous systems.

In this paper we develop and validate a model, based on those employed in Refs. Kremer and Grest (1990); com (a); Baljon et al. (2007), that shows all these characteristics. Much of the physics of amorphous polymer melts and gels is independent of chemical detail Doi and Edwards (1986); Rubinstein and Colby (2003), so properties i) and ii) are captured by using a bead-spring model Kremer and Grest (1990) for the underlying (non-associative) polymer physics. We capture properties (iii)-(v) by employing a hybrid molecular dynamics / Monte Carlo (MD/MC) approach with controlled (specifically, binary) bonding and variable chemical kinetics com (a). The combination of properties (i-v) allows the model to be used to obtain many results unattainable with previous methods. A key advance is that our algorithm is parallelizable. In validating and investigating our model, we separate the effects of thermodynamics, kinetics, and polymer physics on AP network dynamics and mechanics. We show that sticky bonding is mappable to a “mean-field” two-state Arrhenius rate model, but that the kinetic rate constants for SB association and dissociation are affected by the fact that the SMs are embedded in polymers.

One of the key aspects of associating polymer networks, which has rarely been examined (for networks) in the theoretical or simulation literature, is that sticky bond relaxation can be either kinetically limited (i. e. limited by the intrinsic rate of bonding and/or debonding) or diffusion limited (i. e. when the kinetics are so fast that sticky bond scission and recombination events become correlated because newly broken SM pairs tend to recombine before exploring the network cage). Most experimental systems seem to be kinetically limited, and this is the case treated by almost all published analytic theories for reversibly associating networks, including transient network models Tanaka and Edwards (1992); Vaccaro and Marruci (2000). Theories for kinetically limited AP networks Baxandall (1989); Leibler et al. (1991); Rubinstein and Semenov (1998, 2001); Semenov (2006) assume that the sticky bond lifetime $\tau_{sb}$ is so long compared to all “polymeric” relaxation times that it controls all important time scales for network relaxation, and therefore that kinetics only affect key network relaxation times through their effect on $\tau_{sb}$ (or, alternately, as suggested by recent experiments Yount et al. (2005); Loveless et al. (2005), the inverse dissociation rate constant $k_{b}^{-1}$, defined below). We show that the validity of this assumpton is questionable, and that there is a wide parameter space, plausibly accessible in experiments, where it is invalid. Note that published theories de Gennes ; O’Shaughnessy (1987); Friedman and O’Shaughnessy (1988); O’Shaughnessy and Yu (1995); O’Shaughnessy and Vavylonis (1998); Fredrickson and Leibler (1996) for diffusion-limited reactions in polymeric systems, such as those of O’Shaughnessy et. al., either are not immdeiately applicable or have not yet been applied to reversible networks.

Our simulations explore the crossover between the kinetically limited and the diffusion limited cases. We confirm a key prediction of Rubinstein and Semenov Rubinstein and Semenov (1998) on the role of bond recombination on AP dynamics, specifically that recombination effectively renormalizes the SB lifetime in systems where the sticky bonds are sufficiently strong. However, we show that recombination couples interestingly to the diffusive-kinetic crossover in a way not previously predicted.

We extensively examine the dynamics in quiescent, equilibrium systems. One of the more interesting results is that slowing chemical kinetics increases dynamical heterogeneity Kumar and Douglas (2001) in a manner similar to increasing the thermodynamic strength of the sticky bonds. We also examine nonequilibrium ‘chemical dynamics’ (i. e. in systems where the initial sticky bond population is not equilibrated), and nonlinear mechanical properties. The studies of nonequilibrium systems and mechanical properties presented here are limited in number because this paper is intended primarily to illustrate the broad utility of the model.

The organization of the rest of this paper is as follows. In Section II we further motivate, describe and extensively validate our model. We also discuss how it differs from those used in previous simulations of APs. In Section III we present various results for the equilibrium and nonequilibrium physics of thermoreversible AP networks, and compare to theoretical predictions. Finally, Section IV presents a discussion and conclusions. Two Appendices include technical details of the analyses.

II Models and Methods

Of particular interest for the current study, and motivating our modelling approach, are recent experiments performed by the group of Stephen L. Craig Yount et al. (2003, 2005); Loveless et al. (2005); Rotello and S. Thayumanavan (2008); xu ; cra . These have attempted to independently vary thermodynamics and chemical kinetics by making systematic changes in sticky monomer chemistry (based on metal-ligand interactions). Systems with similar static properties show dramatic differences in time-dependent properties that are directly associated with the different kinetics. This effect is quite challenging to capture experimentally because thermodynamics and kinetics are highly correlated for most AP systems (i. e. stronger binding $\leftrightarrow$ slower kinetics Rotello and S. Thayumanavan (2008)), but relatively easy to impose in simulations.

The advantages of including variable kinetics in the model are rather obvious given the above discussion. The advantage of imposing binary bonding (i. e. a SM is at all times bonded to either 0 or 1 other SM) is also relevance to current experiments. Sticky monomers with binary bonding are considered particularly valuable for making thermoplastic elastomers with enhanced melt processability and controllable network architecture Sijbesma et al. (1997); Söntjens et al. (2008). Real examples include multiple hydrogen bonding monomers such as ureidopyrimidinone (UPy), which have highly directional associative interactions, and are of a strength such that the sticky bonds they form constitute a “reversible alternative for the covalent bond” Binder and Zirbs (2007); de Greef and Meijer (2008); Sijbesma et al. (1997); Söntjens et al. (2000, 2008). In the model used here Baljon et al. (2007), sticky bonds differ from covalent bonds only by their reversibility and strength.

II.1 Hybrid MD/MC Simulation Protocol

Our model is built on the framework of the Kremer-Grest bead-spring model Kremer and Grest (1990), which has been extensively validated and is known to capture the key physics of linear homopolymer melts Kremer and Grest (1990); Pütz et al. (2000) as well as permanently crosslinked networks com (c). Each associating polymer chain is linear and contains $N$ beads (monomers). Systems consist of $N_{ch}$ chains, so the total number of monomers is $N_{ch}N$. Periodic boundary conditions are imposed in all three directions, with periods $L_{i}$ along directions $i=x$, $y$, and $z$. Values of $N_{ch}$ range from 700 to 5600; the lowest values are chosen to satisfy $L_{i}>2R_{ee}$ (where $R_{ee}$ is the mean chain end-end distance), preventing self interactions. Remaining leading order finite size effects in networks scale as $(N_{ch}N)^{-1/3}$ Hentschke and Oyen (2005) and should be $\sim 2\%$ for the systems considered here.

All monomers have mass $m$ and interact via the truncated and shifted Lennard-Jones (LJ) potential $U_{LJ}(r)=4u_{0}[(a/r)^{12}-(a/r)^{6}-(a/r_{c})^{12}+(a/r_{c})^{6}]$, where $r_{c}$ is the cutoff radius and $U_{LJ}(r)=0$ for $r>r_{c}$. Covalent bonds between adjacent monomers on a chain are modeled using the finitely extensible nonlinear elastic potential $U_{FENE}(r)=-(1/2)(kR_{0}^{2}){\rm ln}(1-(r/R_{0})^{2})$, with the canonical Kremer and Grest (1990) parameter choices $R_{0}=1.5a$ and $k=30u_{0}/a^{2}$. In this study, following the majority of bead-spring studies on permanently crosslinked systems (e. g. Refs. com (c); Svaneborg et al. (2008)), we employ flexible chains with no angular potential. We express all quantities in units of the LJ bead diameter $a$, intermonomer energy $u_{0}$, and the LJ time $\tau_{LJ}=\sqrt{ma^{2}/u_{0}}$.

All systems have monomer density $\rho=0.85/a^{3}$. We employ two temperatures in this study: $k_{B}T=1.0u_{0}$ and $k_{B}T=0.6u_{0}$. These ambient conditions both correspond to dense polymer melts far above the glass transition temperature $T_{g}$; $T_{g}\simeq 0.35u_{0}/k_{B}$ for $r_{c}=1.5a$ and decreases with decreasing $r_{c}$ Rottler and Robbins (2003); Bennemann et al. (1998). This far above $T_{g}$, melt physics is known to be dominated by the repulsive part of the intermonomer interactions Doi and Edwards (1986); for convenience, and following convention Kremer and Grest (1990), we use purely repulsive LJ interactions with $r_{c}=2^{1/6}a$. However, including attractive interactions by increasing $r_{c}$ is trivial, is important to realistically capture $T$ dependent properties, and will be done in upcoming studies.

All simulations are performed using an enhanced version of the LAMMPS Plimpton (1995) MD code. Newton’s equations of motion are integrated with MD using the velocity Verlet method Frenkel and Smit (2002) and typical timestep $\delta t=.01\tau_{LJ}$ foo (c). A Langevin thermostat Schneider and Stoll (1978) is used to maintain the temperature. The damping time $\tau_{Lang}=10-100\tau_{LJ}$ is larger than the value typically used ($\tau_{Lang}\simeq\tau_{LJ}$) in bead spring studies; this reduces undesirable thermostat-driven effects such as alteration of stress relaxation by suppression of hydrodynamic momentum transfer Pastorino et al. (2007). In this study we employ two “chain lengths”. Most studies are performed at $N=50$, which is at or below best estimates of the entanglement length $50\lesssim N_{e}\lesssim 85$ Everaers et al. (2004); Hoy et al. (2009), so the melts can be fully equilibrated by allowing chains to diffuse several $R_{ee}$ Auhl et al. (2003). $N=50$ is also a convenient choice of chain length because it has been considered in many previous studies. To elucidate the effects of underlying polymeric structure on AP physics, we also consider monomeric melts ($N=1$), that reversibly sticky-bond into dimers.

After the melts are equilibrated, we choose a fraction $c_{st}$ of the monomers to be “sticky”. For the $N=1$ systems these are chosen randomly. For $N=50$ systems, SMs are placed uniformly along chains: at both chain ends and also at internal monomers $iN/(Nc_{st}-1)$, where $i=1,2,...,Nc_{st}-2$. In this study we use $c_{st}=0.08$ (which is comparable to typical experimental values, e. g. Yamauchi et al. (2003); Feldman et al. (2009), so for $N=50$ the SMs are the 1st, 17th, 34th, and 50th monomers in each chain. However, any SM placement can be used, and studies of the effects of altering SM placement at fixed $c_{st}$ are underway; the effects of chemical disorder are known to be significant for stress relaxation Semenov et al. (2007); Feldman et al. (2009).

Sticky monomers are identical to regular monomers, except that they form reversible “sticky” bonds. Figure 1 illustrates the potential energy between sticky monomers as a function of their separation $r$. SMs (like all monomers) always interact via Lennard Jones interactions, whether bonded or not. Bonded SMs additionally interact via the potential $U_{sb}(r,h)$:

[TABLE]

which is based on the standard covalent FENE potential. Here $r_{0}$ represents the equilibrium FENE bond length; $r_{0}\simeq.96a$, i. e. the minimum of the potential $U_{LJ}(r)+U_{FENE}(r)$. The only difference between the sticky and covalent bond potentials is thus an $r$-independent, tunable offset. The same bonding potential was used in Huang et. al.’s studies of equilibrium polymers com (a, d) and a very similar potential was used in Baljon et. al.’s studies of telechelic associating networks Baljon et al. (2007); foo (d). However, our method has several important differences from those of Refs. com (a, d); Baljon et al. (2007) (see Section II.2), so we explain it in detail below.

The potential $U_{sb}(r)$ has several other important features. $h$ represents the sticky binding energy; for $h=0$ a sticky bond can be formed between two monomers separated by $r_{0}$ with no change in energy. The associated force $F_{sb}(r)=-\partial U_{sb}/\partial r$, however, is independent of $h$. Adjusting $h$ is thus a nearly pure way of adjusting the thermodynamics of the sticky bonds without directly altering their “chemical kinetics”, i. e. the rates of formation and dissociation of sticky bonds.

Formation and breaking of sticky bonds is performed using Metropolis Monte Carlo Frenkel and Smit (2002). The change in energy required to form a sticky bond between an unbonded pair of SMs is just $\Delta E(r,h)=U_{sb}(r,h)$, and the energy change to break a sticky bond is $\Delta E(r,h)=-U_{sb}(r,h)$. These are the only MC “moves” used, and the acceptance probability of the moves is set by $\Delta E(r,h)/k_{B}T$. The MC moves are strictly “topological”. All spatial motion of bonded SMs is governed by the sticky bond force $F_{sb}(r)$, along with the other forces from Lennard Jones and covalent FENE interactions, which are all integrated using MD. One potential difficulty is that only bonded SMs “feel” the force $F_{sb}(r)$, so formation/breaking of SMs creates temporal force discontinuities. However, as will be shown below, this does not seem to cause any spurious behavior.

The system sizes and time scales studied here require simulation on parallel computers; 8 to 64 processors are used in a typical simulation. While MD parallelizes very well Plimpton (1995), it has long been noted that MC Esselink et al. (1995) is very difficult to parallelize. We therefore perform hybrid parallel MD / serial MC simulations. MC moves are performed once every $\tau_{0}$ in Lennard Jones time units; $\tau_{0}$ is the MC “timestep”. The MD simulation is paused while the parallel-distributed lists of sticky bonds and SM coordinates are gathered onto one processor. For efficiency, Verlet-style pair neighbor lists of SMs are used and of open SM pairs, only those within $r<R_{0}$ are considered for SB formation. After the SB list is updated, it is distributed back to all processors and the MD simulation resumes. Great care was taken to optimize the MC algorithm to minimize “dead” time on the other processors, but reasonable parallel performance requires $\tau_{0}\gg\delta t$. In this paper, except where otherwise noted, we use $\tau_{0}=\tau_{LJ}$; see also Section II.4.

As mentioned above, current experimental trends favor binary-bonding sticky monomers, We impose binary bonding through a simple restriction on the Monte Carlo routine; sticky bond formation is attempted only for pairs of unbonded SMs. The 1-1 bonding restriction imposed here was also assumed in Ref. Rubinstein and Semenov (1998), which eases comparison of our results to theoretical predictions.

Two further technical details of the MC algorithm are noteworthy: (1) We do not allow any SM pair to both break and form a SB (i. e., to break and recombine) during the same MC step. This is a technical violation of detailed balance, but satisfies the weaker “balance” condition sufficient Manousiouthakis and Deem (1999) for accurate MC simulations. SB recombination is a critical aspect of supramolecular polymer physics Rubinstein and Semenov (1998) and is further discussed in Section III. (2) we do not allow “bond switching” moves within a single MC step. That is, for SMs $V,W,X,Y$, we do not allow moves of the form

[TABLE]

or any other more complicated moves. In addition to being difficult to implement in simulations, such processes are unlikely to occur instantaneously in real systems, in part because of steric constraints. Different but analogous rules, suitably modified to the use of SMs with two binding sites each, were imposed in Refs. com (a, d).

In our model, varying the relative rates of sticky bond formation and dissociation is accomplished by varying $h$. However, the absolute values of the rates depend on a yet unspecified kinetic time scale $\tau_{kin}$. This time can be controlled through the Monte Carlo routine. At each MC timestep (i. e., every $\tau_{0}$), a fraction $f_{MC}$ of unbonded SM pairs (of those within range $r<R_{0}$) and an equal fraction $f_{MC}$ of bonded SM pairs are randomly selected to be considered respectively for sticky bond formation and breaking. We have verified this scheme maintains ‘balance’ Manousiouthakis and Deem (1999) for pairs and triplets of SMs for $.01\leq f_{MC}\leq 1$.

Thus the average time over which each unbonded or bonded SM pair is considered once for (respectively) SB formation or breaking is $\tau_{MC}=f_{MC}^{-1}\tau_{0}$, and the parameter $\tau_{MC}$ effectively controls the “chemical kinetics” of the SBs. For a discussion of why we use $f_{MC}<1$ rather than varying $\tau_{0}$, see Section II.4. Small $\tau_{MC}$ correspond to fast chemical kinetics foo (g), while large $\tau_{MC}$ correspond to slow chemical kinetics.

In Section III.5 we perform mechanical tests on various systems. Two types of tests are perfomed; constant volume deformation and tensile creep. In the constant volume deformation tests, $L_{z}$ is increased at a true strain rate $\dot{\epsilon}=\dot{L_{z}}/L_{z}$, and $L_{x}$ and $L_{y}$ are adjusted to maintain constant volume. In the creep tests, a constant (small) stress difference $\sigma_{creep}$ (relative to the equilibrium hydrostatic pressure in the quiescent state, which is positive for repulsive LJ interactions) is applied along the $z$ direction using a Nose Hoover barostat Frenkel and Smit (2002). This smaller $|\sigma_{z}|$ produces tensile creep. Both types of tests use $\tau_{0}=.2\tau_{LJ}$ to minimize systematic errors.

II.2 Comparison to Previous Simulation Protocols

It is worthwhile to compare the simulation method and ambient conditions described above in the context of previous AP simulation studies. The use of a hybrid MD/MC method is a powerful advantage. Pure Monte Carlo (MC) simulations have been performed with lattice Nguyen-Misra and Mattice (1995); Y.Liu and Pandey (1996); Kumar and Panagiotopoulos (1999); Kumar and Douglas (2001); Wang et al. (2008) and off-lattice Groot and Agterof (1994); Loverde et al. (2005) models. These are very effective at studying static properties like percolative gelation and (in the case of solutions) phase separation, but have limited ability to capture the complex, collective dynamics which occur in bulk polymers, and thus lack properties (i) and (ii). For example, MC can not, even in principle Milchev et al. (1999), capture hydrodynamic effects, which are expected to play an important role whenever momentum transfer is important (e. g. in relaxation of highly stressed systems).

Pure molecular dynamics (MD) studies Khalatur et al. (1999); Ayyagari et al. (2004); Guo and Luijten (2005); Manassero et al. (2005); com (b); Goswami et al. (2008) have been used to study static and dynamic properties. While better able to capture dynamics and nonequilibrium phenomena than MC, MD studies can not naturally implement controllable sticky bonding topology. Also, MD studies cannot easily impose any control of chemical kinetics without resorting to costly, chemically realistic models. For example, Padding and Boek Padding and Boek (2004) studied systems intermediate between ours and those studied by Huang et. al.; a fraction $c_{st}<1$ of their monomers were allowed to form linear equilibrium poiymers, but the FENE-C sticky bonding potential Kröger (2004) used did not allow for variable kinetics. Thus, in practice, typical MD studies lack properties (iii) and (iv).

The previous works most closely related to the present method are Refs. Baljon et al. (2007); com (a, d), who also used hybrid MC/MD with the same $U_{sb}(r)$ (Eq. 1). Huang et. al. com (a, d) also used variable kinetics, but studied equilbrium linear polymers with $c_{st}=1$ rather than network-forming APs with $c_{st}\ll 1$. Details of the Huang et. al. method are discussed extensively in Ref. com (a). The key differences of our method from Ref. Baljon et al. (2007) are the imposition of binary bonding and the use of variable $\tau_{MC}$ (they used only one $\tau_{MC}=0.2\tau_{LJ}$). Another difference was that Refs. Baljon et al. (2007); com (a, d) all used a much stronger thermostat, giving overdamped (Brownian) dynamics.

Many previous studies have used nonspecific (e. g. strengthened attractive Lennard-Jones or Coulombic) interactions which allow SMs to form arbitarily many simultaneous SBs Nguyen-Misra and Mattice (1995); Khalatur et al. (1999); Kumar and Panagiotopoulos (1999); Kumar and Douglas (2001); Ayyagari et al. (2004); Guo and Luijten (2005); com (b); Manassero et al. (2005); Goswami et al. (2008). This results in formation of interesting structures such as micelles and micelle-bridge networks, which occur in real AP systems such as associating ionomers (see e. g. Ref. com (e)). In contrast, the (experimental) APs we wish to model tend to form networks more like classical rubbers.

Most previous studies Groot and Agterof (1994); Y.Liu and Pandey (1996); Khalatur et al. (1999); Kumar and Douglas (2001); Loverde et al. (2005); Baljon et al. (2007) varied temperature $T$ at fixed SM bonding strength. This does not isolate the effects of $T$ on associative bonding from its other effects such as the dynamical slowdown which occurs in normal (non-associative) polymers. To get a full picture of AP physics, one should vary both $h$ and $T$ independently com (a). We follow this approach.

Other differences from previous simulation studies are more associated with the systems employed than the methods applied. Many studies have considered only telechelic chains Khalatur et al. (1999); com (b); Loverde et al. (2005); Guo and Luijten (2005); Ayyagari et al. (2004); Manassero et al. (2005); Baljon et al. (2007); Goswami et al. (2008). Telechelics are appealing in their simplicity, but their network-forming abilities are naturally limited; for binary bonding, at least 3 SMs/chain are required to form good networks. Weakly entangled chains ($N\sim N_{e}$) may be ideal de Greef and Meijer (2008) for technological goals such as enhanced melt processability at high $T$ and network strength at low $T$. The majority of previous studies have employed extremely short $N\ll N_{e}$ chains Manassero et al. (2005); Y.Liu and Pandey (1996); Khalatur et al. (1999); Ayyagari et al. (2004); Guo and Luijten (2005); Baljon et al. (2007), but we consider systems with $N\sim N_{e}$. Finally, the majority of previous studies have focused on small $\rho$ corresponding to solutions Nguyen-Misra and Mattice (1995); Khalatur et al. (1999); Kumar and Panagiotopoulos (1999); Kumar and Douglas (2001); Guo and Luijten (2005); Baljon et al. (2007); com (a, d). AP solutions exhibit a wide range of intriguing phenomena, in particular competition between gelation and phase separation Groot and Agterof (1994); Y.Liu and Pandey (1996), which, however, we do not wish to consider here. In addition, the presence of solvent can dramatically weaken the effective strength of sticky bonds in real systems Söntjens et al. (2000); Yamauchi et al. (2003); de Greef and Meijer (2008); this effect is beyond the scope of our model. We therefore focus on systems with $\rho$ corresponding to a dense pure melt with no solvent.

Mappings of the bead-spring model to real, dense polymer melts Kremer and Grest (1990) produce different $\tau_{LJ}$ in the range $10^{-10.5\pm 1.5}s$. Present day computers can achieve runs (for the system sizes used here) of up to $\sim 10^{7}\tau_{LJ}\sim 10^{-3.5\pm 1.5}s$, but runs this long can not be performed over a broad parameter space. In contrast, sticky bond lifetimes in experimental systems are typically at least $10^{-4}s$, and often many orders of magnitude longer Sijbesma et al. (1997); Söntjens et al. (2000); Yount et al. (2005). Thus any attempt to capture specific SM chemistries and at the same time use systems large enough to study bulk dynamics and mechanical properties would exceed the capabilities of present day supercomputers foo (e). Coarse-grained modelling with the goal of studying the dynamics of AP systems by analogy is the only currently feasible approach for bulk systems, so we make no attempt to mimic specific chemistries. The only published simulations of which we are aware that model AP networks with specific chemistries Wang et al. (2008); Genesky et al. (2008) are pure Monte Carlo studies that used a very coarse-grained (lattice) bond-fluctuation model Carmesin and Kremer (1988) and focused on static properties.

II.3 Static Properties: Validation of Hybrid MD/MC Method

As discussed above, systems contain a total of $N_{st}=N_{ch}Nc_{st}$ sticky monomers. Due to the binary bonding rules, the maximum number of sticky bonds that can exist in the system at any given time is $N_{st}/2$. If the probability that an SM is bound into an SB is $p_{active}$, then the total number of SBs in the system is $N_{st}p_{active}/2$. If $A$ represents an unbound SM and $A_{2}$ represents a bound SM pair, these factors define the concentrations

[TABLE]

where square brackets denote concentrations. If the equilibrium value of $p_{active}$ is $p^{*}$, then the equilibrium constant for SB association is defined (by the law of mass action for the reaction $A+A\leftrightarrow A_{2}$) as

[TABLE]

for binary bonding.

Figure 2 shows simulation data in which $p^{*}$ was evaluated from equilibrated simulations at fixed $h$ and $K_{eq}$ obtained from Eq. 4. Circles show values of $K_{eq}$ for $N=1$ and $N=50$ systems. As expected, $K_{eq}\sim exp(h/k_{B}T)$. The data shown are for $\tau_{MC}=1.0\tau_{LJ}$, but we have verified that $p^{*}$ is independent of $\tau_{MC}$ (to within statistical errors) for all $h$ tested, over the range $\tau_{LJ}\leq\tau_{MC}\leq 100\tau_{LJ}$. Because there is an entropy cost $\sim k_{B}T$ to form a SB, few SBs form for $h<2u_{0}$. As $h/u_{0}$ ranges from $2$ to $17.5$, the equilibrium constant $K_{eq}$ varies over more than six orders of magnitude, from $0.96$ to $3.9\cdot 10^{6}$. This is a wider range of $h$ and $K_{eq}$ than considered in previous simulation studies. A standard Pandey and Landau (1985) finite size scaling analysis of the percolation gel transition is given in Appendix A. For $k_{B}T=1.0u_{0}$, percolation occurs at $h=h_{perc}=4.25u_{0}$, so we consider values of $h$ up to $\sim 4$ times above the gelation transition.

Note that the $\tau_{MC}$-independence of $p^{*}$ allows systems to be equilibrated efficiently using a low $\tau_{MC}=\tau_{LJ}$. Higher values of $h$ (for polymeric systems) are impossible to equilibrate on present-day computers with our current method; equilibration is discussed further in Section III. However, the highest values of $K_{eq}$ considered here are comparable to those observed in some experiments on multiple-H-bonding SMs Söntjens et al. (2000); Rotello and S. Thayumanavan (2008).

Data from multiple system sizes are also useful in further validating the simulation model. Ben-Naim and Krapivsky have pointed out that systems which reversibly polymerize undergo a nonthermodynamic gelation transition Ben-Naim and Krapivsky (2008) when the fragmentation (in our case, SB breaking) process is too weak. The average number of clusters (aggegrates) at any given time is $N_{agg}\equiv N_{ch}/N_{n}$, where $N_{n}$ is the number-averaged cluster size (Appendix A). $N_{agg}\equiv N_{ch}$ in the absence of sticky bonding and $N_{agg}\to 1$ in the limit of large $h$, because all the chains combine into a single network (as in an ideal rubber). Our systems, in the terms of Ref. Ben-Naim and Krapivsky (2008), are “thermodynamic” if and only if: (1) $N_{agg}$ is linearly proportional to $N_{ch}$ below percolation (i. e. for $h<h_{perc}$) and (2) the probability distibution of cluster sizes $P(M)$ (Appendix A) is independent of $\tau_{MC}$. An arbitrary simulation method will not necessarily display a ‘thermodynamic’ gel transition; failure to do this would be a serious flaw according to our goals. We therefore have verified that our model satisfies conditions (1) and (2) for $\tau_{LJ}\leq\tau_{MC}\leq 100\tau_{LJ}$, and therefore properly captures reversible gelation. Satisfaction of these conditions appears equivalent to the above-verified condition that $p^{*}$ is independent of $\tau_{MC}$ foo (g).

Figure 3 shows data for $P(M)$ at $h=4u_{0}$ and $k_{B}T=1.0u_{0}$ (i. e. just below percolation). The collapse of the data shows Ben-Naim and Krapivsky (2008) that cluster formation/dissociation is an equilibrium processes and supports our arguments that the algorithm satisfies detailed balance for the range of $\tau_{MC}$ considered here. Also, $P(M)$ shows some interesting properties which demonstrate that our modelled systems form good (rubber-like) networks. The line shows a fit to a $P\propto(M)^{-5/4}$ power law, which is consistent with the fractal dimension $D_{frac}=4$ of aggregrates and the expected power law $ln(P)\sim-(1+1/D_{frac})ln(M)$ for networks Milchev et al. (2000). In contrast, dense telechelic systems have an exponential $P(M)\sim exp(-M/<M>)$ distribution. The absence of any large exponential contribution in our $P(M)$ at large $M$ indicates that long linear clusters are not common. Therefore, though our parent chains only contain 4 SMs each, we are confident that that is enough to accurately capture AP network physics.

II.4 SB Dynamics and Two-State Model

Figure 4 shows simulation results for the average sticky bond lifetime, $\tau_{sb}$, in quiescent systems at chemical equilibrium. Simple thermal activation of SB dissociation would suggest exponential behavior, $\tau_{sb}^{-1}\propto exp(-h/k_{B}T)$. In fact the results are markedly nonexponential. Interestingly, SB lifetimes in polymeric systems are (apparently) always lower than those in dimer-forming systems. This is consistent with differences in chain connectivity; SBs embedded in polymers experience additional ‘pulling’ forces due to transmission of the random thermal forces (which produce diffusive motion) through covalent bonds along their parent chains. Additional reductions in $\tau_{sb}$ could potentially arise from increased steric hindrance to bonding for embedded SMs. While this “polymeric” effect on $\tau_{sb}$ should dependent sensitively on $N$, $c_{st}$ and $T$, to our knowledge it is not included in any theories for AP networks.

The simulation data in Figures 2 and 4 can be better understood by mapping the Monte Carlo procedure and $U_{sb}(r,h)$ onto a two state Arrhenius model for sticky bonding. The model is depicted in Figure 5. Bonded SM pairs are assumed to have an energy $-h$, unbonded SMs have zero energy, and we introduce an $h$-dependent barrier $\delta(h)$.

The Monte Carlo rules described in Section II.1 allow us to assume that sticky bond formation obeys second order kinetics and dissociation obeys first order chemical kinetics, as they should as long as $\rho c_{st}\ll 1/a^{3}$ Atkins (1978); Rubinstein and Semenov (1998). The SB formation/dissociation process can be represented as the chemical reaction

[TABLE]

where $k_{f}$ and $k_{b}$ are the rate constants for SB formation and dissociation. Then the equation for chemical equilibrium is

[TABLE]

In the Arrhenius two state model the rate constants are given by:

[TABLE]

where $\alpha$ and $\beta$ are constants with dimensions of volume$\times$frequency and frequency, respectively. Note that the above is a “mean field” model Cates (1987) in that it ignores correlations between sticky monomers (i. e. concentration fluctuations). Thus $k_{f}$ and $k_{b}$ (and especially $\alpha$ and $\beta$) will in general depend Loverde et al. (2005) on $N$, $\rho$, $c_{st}$, and (through second order effects such as the variation of $\rho$ at fixed pressure) $T$.

In thermal equilbrium, Eq. 6 gives the equilibrium constant

[TABLE]

Eq. 8 fits simulation results for $K_{eq}$ very well, as shown in Fig. 2. In AP networks at even higher values of $h$, Eq. 8 should fail due to ‘trapped’ open SMs Loverde et al. (2005) that cannot find partners, but this effect is negligible for the systems considered here.

We now compare two state model predictions to simulation data and map the latter to the former. In the two state model, the mean SB lifetime $\tau_{sb}$ is just $\tau_{sb}=k_{b}^{-1}$. Similarly, the probability that an unbonded pair in the ‘2A’ state (Fig. 5) will jump over the barrier is just $exp(-\delta(h)/k_{B}T)$. This is also the success rate $S_{MC}\equiv k_{f}/\alpha$ for Monte Carlo SB formation attempts, so $\delta(h)$ can be directly measured from the simulations: $\delta(h)/k_{B}T=-ln(S_{MC})$.

In Figure 6, panel (a) shows simulation results for $\delta(h)$ and panel (b) shows simulation results for $\tau_{sb}^{-1}\exp(\delta(h)/k_{B}T)$. The latter shows that the perfect exponential decay expected from Eq. 7, $\tau_{sb}^{-1}\exp{(\delta(h)/k_{B}T)}=k_{b}=\beta\exp{(-h/k_{B}T)}$, is actually observed. Note that this Arrhenius behavior was in no way imposed; it emerges naturally, showing the utility of the two state model in understanding the behavior of our simulations.

The parameters $\alpha$ and $\beta$ can be extracted from the data in Figs. 2 and 6. For $\tau_{MC}=\tau_{LJ}$, $\alpha=6.3a^{3}/\tau_{LJ}$ and $\beta=24/\tau_{LJ}$ for dimers, while $\alpha=4.0a^{3}/\tau_{LJ}$ and $\beta=41/\tau_{LJ}$ for $N=50$ chains. The smaller $\alpha$ measured for polymer-embedded SMs is consistent with the above-hypothesized increased steric constraints. The large rate constant $\beta=41/\tau_{LJ}$ indicates a potential problem with the simulations. $\beta$ is an effective “attempt frequency” for breaking sticky bonds, which implies that the MC timestep $\tau_{0}$ should be small compared to $\beta^{-1}$. Larger $\tau_{0}$ will in principle produce systematic errors. The data shown above are for $\tau_{0}=1.0\tau_{LJ}$, which is large compared to $\beta^{-1}$. Simply reducing $\tau_{0}$ is problematic because it sharply reduces the parallel efficiency of the simulations.

However, we have used values of $\tau_{0}$ as small as $.05\tau_{LJ}$, and find that all errors produced by using $\tau_{0}=1.0\tau_{LJ}$ are small in quiescent systems at equilibrium; for example, the systematic error in $\tau_{sb}$ at $h=10u_{0}$ is about 1%. While the errors in dynamical properties are somewhat larger at small $h$ ($h\ll 10u_{0}$), in this paper we focus on dynamics for $h\geq 10u_{0}$ and use $\tau_{0}=1.0\tau_{LJ}$. In all cases, all differences produced by smaller $\tau_{0}$ are small compared to the differences between systems contrasted in Section III, and comparable to our statistical errors, i. e. $\sim 1\%$. For nonequilibrium systems, however, systematic errors are larger. Thus all nonequilibrium and mechanical-property tests in this paper are performed using $\tau_{0}\leq 0.2\tau_{LJ}$.

In summary, to within our noise, increasing $\tau_{MC}$ leaves the static properties of our model AP networks (Figure 2) unchanged, and changes the sticky bonding dynamics (Figure 6a-b) only through the prefactor $\tau_{sb}\propto\tau_{MC}^{-1}$. The role of $\tau_{MC}$ in the dynamics therefore appears in the rate constants $\alpha$ and $\beta$, which are also proportional to $\tau_{MC}^{-1}$. Increasing $\tau_{MC}$ slows down the chemical kinetics of the SBs (both formation and dissociation) relative to the underlying polymeric time scales, while leaving the thermodynamics unchanged. This is why we claim our model can separate thermodynamics and kinetics. The variation of $h$ and $\tau_{MC}$ employed here may be thought of as corresponding to “scanning” across chemically different sticky monomers. Given the time scale problem mentioned above, this scanning is only qualitative. However, we show below that it is very useful in understanding AP systems.

III Results

Previous work has shown Kumar and Douglas (2001); Baljon et al. (2007) that the most dramatic changes in dynamics, our primary interest, take place not at $h_{perc}$ but rather at considerably higher $h$. The rest of this paper considers systems with $h\gg h_{perc}$ and $p^{*}\gtrsim 0.95$. This is the “physical gel” regime Kumar and Panagiotopoulos (1999) where nearly all chains are (at any moment) part of a single aggregate. A “snapshot” of a physical gel looks much like a crosslinked rubber, yet chains are delocalized and the system can flow at long times. One of the most interesting properties of physical gels is their transition to chemical gels as $h$ increases or $T$ decreases. In this “physical-chemical gel transition” (PCGT), chains become localized Kumar and Douglas (2001); Baljon et al. (2007) in a manner analogous to the “caging” effect produced upon cooling fragile glass-forming systems Donati et al. (1999); foo (f).

For the $N$, $T$, and $c_{st}$ considered here, the PCGT occurs foo (g) at bonding strength $h_{PCGT}>17.5u_{0}$. The broad range ($5u_{0}\lesssim h\lesssim 17u_{0}$) between the percolation (Appendix A) and localization transitions is consistent with the findings of Kumar and Douglas Kumar and Douglas (2001) as well as Baljon et. al. Baljon et al. (2007), who both, however, used constant SB strength and varied $T$. The broad range is not dependent on having only a few sticky monomers per chain, although increasing $Nc_{st}$ at fixed $\rho$ will broaden the range by lowering $h_{perc}$. Here we focus on values of $h$ which are well below $h_{PCGT}$, and thus “in the middle” of the physical gel regime.

Another of the key features of physical AP gels is sticky bond recombination. The concentration $\rho c_{st}(1-p^{*})$ of free SMs is small. Moreover, the motion of free SMs is constrained by their (transiently but usually) bonded intrachain neighbors. Thus SM pairs tend to recombine after SB-dissociation events. This leads to to a second characteristic timescale for individual sticky bonds; in addition to the “bare” lifetime $\tau_{sb}$, there is Rubinstein and Semenov (1998) a larger, “effective” SB lifetime $\tau^{*}$, which can be thought of as the average time for initially bonded SMs to “separate” (i. e. no longer recombine) as opposed to merely debond. It is of interest because rheological experiments typically measure $\tau^{*}$; $\tau_{sb}$ is more difficult to access Cates (1987); Rotello and S. Thayumanavan (2008); Feldman et al. (2009). Values for $\tau_{sb}$ and $\tau^{*}$ (defined more specifically in Appendix B and discussed further in Section III.2) for a wide variety of systems are given in Table 1. We have already shown how $\tau_{sb}$ is affected by polymer physics - indirectly through covalent backbone bonds. Now we study the ways in which SB recombination influences and is influenced by the interplay of SB thermodynamics, SB kinetics, and polymer physics. We perform our study in terms of measurements of diffusion, $\tau^{*}$, dynamical heterogeneity, nonequilibrium chemical dynamics, and nonlinear mechanical properties. All results presented below are for systems that were first equilibrated for many $\tau^{*}$. As will be shown, these are best understood by determining whether SB recombination is diffusion-limited or kinetically limited.

III.1 Diffusion

The effect of varying different thermodynamic and kinetic parameters on monomer diffusion (mean squared displacement $<(\delta\vec{r})^{2}(t)>$) is shown in Figure 7. Panel (a) shows the variation as $h$ is increased at $\tau_{MC}=1.0\tau_{LJ}$ and $k_{B}T=1.0u_{0}$. At short times ($t\ll\tau_{sb}$), results for different values of $h$ collapse, showing (as expected) that sticky bonding has little effect on diffusion on these time scales. At larger times ($t\gtrsim\tau_{sb}$) results show a progressive localization and ‘caging’ effect, similar to that described in Refs. Kumar and Douglas (2001); Baljon et al. (2007), as sticky bond strength increases. At $h=10u_{0}$, little localization occurs because Rubinstein and Semenov (1998) $\tau_{sb}$ is less than the Rouse time of the chains in the absence of sticky bonding ($\tau_{R}\simeq 2.6\cdot 10^{3}\tau_{LJ}$). As $\tau_{sb}$ increases with increasing $h$, the curves develop a “shoulder” which illustrate the temporary caging associated with physical gels. This temporary cage becomes permanent as $\tau_{sb}\to\infty$ (as in a classical crosslinked rubber). Data in panel (a) support our earlier statement that this occurs for some $h\gtrsim 17u_{0}$.

Panels (b-c) show a pair of interesting effects. First, for $h=11.25u_{0}$ and $k_{B}T=1.0u_{0}$, increasing $\tau_{MC}$ has the same qualitative effect as increasing $h$ at fixed $\tau_{MC}$. Data for $<(\delta\vec{r})^{2}(t)>$ collapse for $t$ less than the smallest $\tau_{sb}$ (i. e. $\tau_{sb}$ for the lowest $\tau_{MC}$). For longer times, the data develops a shoulder which increases in width as $\tau_{MC}$ increases. Data for an equilibrated system with MC deactivated (i. e. $\tau_{MC}=\infty$) shows “chemical gel” (ideal rubber) behavior; chains are permanently localized.

Second, data from systems with the same “SB thermodynamics” (i. e. the ratio $h/k_{B}T=11.25$) but different ambient conditions ($h=6.75u_{0}$ and $k_{B}T=0.6u_{0}$) shows interesting contrasts which illustrate the interplay of SB dissociation and underlying polymer physics. For $t\lesssim\tau_{sb}$, data for $<(\delta\vec{r})^{2}(t)>$ still collapse, but data from the lower-$T$ systems collapse on a lower value. This is not at all surprising, as polymeric diffusion is well known to slow with decreasing $T$. However, though $h/k_{B}T$ is the same, values of $\tau_{sb}$ are smaller for the lower-$T$ systems, perhaps because $h$ is smaller and SB breaking is favorable at smaller $r$ (see Fig. 1) foo (h). Thus the diffusion data actually cross over at intermediate time scales (panel c) and the lower $T$ systems show greater mobility at fixed $h/k_{B}T$, a most unusual state of affairs. While the case presented here is somewhat artificial because in a real polymer melt $\rho$ would decrease with $T$ and lead to further diffusive slowdown, we believe the point that varying $T$ at fixed SB thermodynamics should change relaxation at different timescales differently should be generally valid. For example, the frequency ($\omega$) dependence of the dynamical moduli $G(\omega;T)$ Ferry (1980) should change with $T$ in nontrivial ways. In other words, time-temperature superposition should be violated.

It is useful to relate the mean squared displacement to the cage size $a_{cage}$ and “escape parameter” $f(t)$ using the definition

[TABLE]

In the $\tau_{sb}\to\infty$ limit, $a_{cage}^{3}\sim(\rho c_{st})^{-1}$ is the volume explored by sticky monomers foo (i). The “chemical gel” time $t_{chem}$ is the time at which $f(t)$ approaches unity in this limit. Data for the $\tau_{MC}=\infty$ system in Figure 7(b) (with $h=11.25u_{0}$, $N=50$, and $c_{st}=.08$) shows that $a_{cage}^{2}\simeq 23a^{2}$ and $t_{chem}\sim 10^{4.5}\tau_{LJ}$. For finite $\tau_{sb}$, one can define $t_{cage}$ as a “caging” time describing the (de)localization of SMs Kumar and Douglas (2001); Baljon et al. (2007); $f(t)$ then has the general form $f(0)=0$, $f(t)\sim 1$ for $t\sim t_{cage}$, and $f(t)>1$ for $t>t_{cage}$.

The monomeric diffusion constant $D$, as measured by $lim_{t\to\infty}<(\delta\vec{r})^{2}(t)>\sim 6Dt$, should vary inversely with some “long” characteristic time $\tau_{long}$ of the system, roughly defined as the time for chains to diffuse by their end-end distance. Candidates for $\tau_{long}$ include $\tau_{sb}$ and $\tau^{*}$. Ref. Rubinstein and Semenov (1998) predicts $\tau_{long}\propto\tau_{sb}$ for weakly binding physical gels and $\tau_{long}\propto\tau^{*}$ in the strong-binding (near-chemical) limit. Figure 8 shows results for $D\tau_{sb}$, $D\tau^{star}$, and $8\cdot 10^{4}\tau_{LJ}D$ from Table 1 for $h=11.25u_{0}$, $k_{B}T=1.0u_{0}$ systems, over a wide range of $\tau_{MC}$. $D$ decreases with increasing $\tau_{MC}$ slower than both $(\tau^{*})^{-1}$ and $\tau_{sb}^{-1}$, but it tracks the former more closely than the latter. Thus results for this value of $h$ are apparently intermediate between the “weak” and “strong” physical gel limits described in Ref. Rubinstein and Semenov (1998).

III.2 Crossover between Diffusion-Limited and Kinetially Limited SB recombination

Table 1 shows $\tau^{*}$ and $\tau^{*}/\tau_{sb}$ for all investigated systems. As expected, increasing $h$ at fixed $\tau_{MC}$ increases $\tau^{*}/\tau_{sb}$, because for fixed kinetics recombination is more likely for thermodynamically stronger SBs. There are several possible regimes of possible relations between $\tau^{*}$ and $\tau_{sb}$ that can be related to diffusion (specifically, $<(\delta\vec{r})^{2}(t)>$) on intermediate timescales. Systems with $\tau_{sb}\gg min(t_{cage},t_{chem})$ will exhibit kinetically limited sticky bond recombination (KL); in this regime $\tau^{*}/\tau_{sb}$ is predicted to be constant Rubinstein and Semenov (1998). However, if $\tau_{sb}\ll\tau_{cage}$, recombination will be dominated by “correlated” recombinations of SM pairs that have recently dissociated and have not had time to fully diffuse away from one another. This is diffusion limited sticky bond recombination (DL). To our knowledge, the DL regime and especially the crossover between DL and KL have not been previously studied for AP networks.

Figure 9 shows the variation of $\tau^{*}/\tau_{sb}$ with chemical kinetics for $h=11.25u_{0}$ systems at $k_{B}T=1.0u_{0}$. Over a range of two orders of magnitude in $\tau_{MC}$, the data are well fit by the equation

[TABLE]

where $C$ is the probability of recombination in the KL limit and $K$ is the contribution from diffusion-limited recombination. The exact form of Eq. 10 is not of great consequence; what matters is the broad crossover between regimes and the large change of $\tau^{*}/\tau_{sb}$ as a function of kinetics. Nevertheless, since $\tau_{sb}=k_{b}^{-1}$, Eq. 10 can be interestingly rewritten

[TABLE]

The significance of Eq. 11 is its prediction of a nonlinear dependence of $k_{b}\tau^{*}$ on the rate constant $k_{b}$ for dissociation; $\tau_{MC}^{-x}\propto k_{b}^{x}$ (see Section II.4).

$C$ and $K$ are of course not universal constants, but will depend on $h$, $\rho$, $c_{st}$, $T$, and $N$. In practice, one would expect $K\tau_{sb}^{x}\ll C$ when $\tau_{sb}\gg min(t_{cage},t_{chem})$. More physically, the condition $K\tau_{sb}^{x}\ll C$ defines the KL regime, where sticky bond reactions become “mean-field” in the sense of Cates Cates (1987). The data in Figure 9 show that one can (at least in our model systems) move from the KL to the DL regimes simply by speeding up the chemical kinetics, if one is in the regime where $\tau_{sb}$ is comparable to the underlying polymer relaxation times such as $t_{chem}$.

For our systems, values of $C$ are close to values of $\tau^{*}/\tau_{sb}$ in our “kinetically slow” ($\tau_{MC}=100\tau_{LJ}$) systems, cf. Table 1. $C$ increases with $h$ (qualitatively) as predicted by Rubinstein and Semenov Rubinstein and Semenov (1998). However, while our kinetically slow systems all have $C<2$, the “strong physical gel” theory in Ref. Rubinstein and Semenov (1998) assumes $C\gg 1$, so we defer a detailed comparison to that theory to later work. Here we merely make the positive observation that the basic prediction Rubinstein and Semenov (1998) of effective SB lifetime renormalization ($\tau_{sb}\to\tau^{*}$) works well at these relatively small $C$ and (somewhat surprisingly) over the entire studied KL$\to$DL crossover regime. The renormalization $\tau_{sb}\to\tau^{*}$ accurately captures effective SB dissociation over a very broad parameter space foo (i).

On the other hand, an observation apparent from Fig. 9 is that SB recombination is in general only partly ‘polymeric’ in nature. In the limit of fast kinetics, values for $\tau^{*}/\tau$ in $N=1$ systems are nearly as high as for $N=50$ systems. However, the lack of a network prevents any caging effects, and so $\tau^{*}/\tau_{sb}$ decreases much faster with increasing $\tau_{MC}$, approaching unity (regardless of the value of $h$) at $\tau_{MC}=100\tau_{LJ}$. This illustrates that the crossover from DL to KL is analogous to a crossover between ‘dimeric’ and ‘polymeric’ recombination; in other words, chain connectivity (i. e. covalent bonding) becomes increasingly important as kinetics are slowed. While the ‘dimeric’ contribution to SB recombination cannot be simply “subtracted out” due to the different $\delta(h)$, these ‘dimeric-vs-polymeric’ effects on SB recombination have been neglected by previous theories foo (j).

III.3 Recombination and Dynamical Heterogeneity

Figure 10a shows simulation results for the SB recombination probability $P_{recomb}(t)$ (defined in Appendix B) for $h=11.25u_{0}$ systems with different chemical kinetics. Although there is a small peak in $P_{recomb}$ (not displayed and low compared to the peaks shown in Fig. 10a) at very small times $t\sim\tau_{MC}$, the $t^{-5/4}$ behavior predicted O’Shaughnessy and Yu (1995) for the extreme diffusion-limited case is not found. This indicates none of our systems have “too-fast” kinetics com (a). For $t\gg\tau_{0}$, our results have the interesting form $P_{recomb}(t)\simeq\exp(-t/\tau^{*})-\exp(-t/\tau_{sb})$; a comparison to actual data for $\tau_{MC}=10\tau_{LJ}$ is shown.

This form of $P_{recomb}(t)$ has a maximum at the “delocalization” time

[TABLE]

where $y=\tau^{*}/\tau_{sb}$. Bonded SM pairs trend towards moving away from each other after $\tau_{deloc}$. Furthermore,

[TABLE]

which becomes small for $y\gg 1$. This suggests $\tau_{deloc}$ (in addition to $\tau^{*}$) might be a key relaxation time in systems with large $y$. However, this is speculative and needs further verification.

A useful measure of relaxation in complex fluids is the “non-Gaussian” parameter

[TABLE]

which is zero for normal diffusion and positive for systems where some particles move anomalously fast Guo and Luijten (2005), particularly for “hopping” type motion. $\Lambda$ has been shown to be relevant to the structural relaxation of supercooled liquids and dynamical heterogenity Donati et al. (1999). The time $t_{deloc}$ at which $\Lambda$ is maximized and the maximum value $\Lambda_{max}=\Lambda(t_{deloc})$ both increase with decreasing $T$ in various systems, including associating polymers Kumar and Douglas (2001); Guo and Luijten (2005); Baljon et al. (2007), as localization increases. $t_{deloc}$ may be regarded as a crossover time after which the system begins to show liquidlike behavior.

Figure 10b shows the effect of kinetics on $\Lambda(t)$. The effect of slowing kinetics at fixed $h/k_{B}T$ is similar to the effect of increasing $h/k_{B}T$ observed in previous studies Kumar and Douglas (2001); Baljon et al. (2007). It is interesting that increasing $\tau_{MC}$ increases dynamical heterogeneity. The probable reason is that increasing $\tau_{MC}$, even though it leaves $p^{*}$ unaffected, decreases the likelihood of multiple closed SBs on the same chain breaking within a short time period. This is consistent with the idea Leibler et al. (1991) that coherent breaking of nearby SBs along a chain eases large-scale motion. The increasing dynamical heterogeneity with increasing $t_{cage}$ is consistent with other results showing $\Lambda_{max}$ increases as localization “transitions” are approached, e. g. stretched-exponential relaxation of finite clusters Ayyagari et al. (2004).

The data in Figs. 9-10 also clearly show that $t_{deloc}\simeq\tau_{deloc}$ in systems where recombination is likely (i. e. when $\tau^{*}/\tau_{sb}$ is large compared to 1), and thus that delocalization is closely related to individual sticky bonds finding new partners in a “hopping” type motion. However, these delocalization times are large compared to $t_{cage}$. This is not surprising, as full delocalization should occur only when chains have lost all memory of their initial SB topology; this “memory” time is inherently polymeric in that it must increase with increasing $Nc_{st}$, similarly to a Rouse or reptation time Doi and Edwards (1986). Interestingly, the peaks of $\Lambda$ are broader than those of $P_{recomb}$. This also likely arises either from cluster effects Rubinstein and Semenov (1998); Ayyagari et al. (2004) or other underlying many-SM phenomena that ultimately arise from the ‘polymer physics’, i. e. the covalent connectivity of the parent chains.

III.4 Nonequilibrium Chemical Dynamics

An important feature of our model is its ability to accurately capture the dynamics of systems in which the sticky bonds are not in thermal equilibrium. The evolution of SB concentration is, following Equation 5, given by

[TABLE]

which after plugging into Eq. 3 and simplifying becomes

[TABLE]

Equation 16 has an analytic solution. For the special initial condition $p_{active}(t)=0$ at $t=0$, the solution is

[TABLE]

where $z=\rho c_{st}k_{f}$ and $d=1+k_{b}/4z$.

Figure 11 compares this analytic prediction to simulation results for $p_{active}(t)$ upon activation of sticky bonding for two systems with the same value of $h/k_{B}T$ but different values of $h$ and $k_{B}T$. Values of $z$ and $d$ in Eq. 17 are taken from fit values of $\alpha$, $\beta$ and the measured value of $\delta(h)$ as reported in Section II; note that these vary somewhat with $T$, giving different $p^{*}$ at the same $h/k_{B}T$. Data agree excellently with predictions at short and long times. The merely qualitative agreement at intermediate times is no cause for concern, but is an interesting ‘feature’, because Eq. 16 ignores all physics arising from the important fact that the sticky monomers are embedded, at a concentration $c_{st}$, in chains of length $N$, in a dense polymer melt. The slower convergence of simulation results for $p_{active}(t)$ relative to the prediction of Eq. 17 is consistent with such polymeric effects; better agreement is observed for dimer systems. As expected, the polymeric slowdown is greater at lower $T$.

The results above demonstrate the ability of our method to capture the effect of polymer physics on nonequilibrium “chemical dynamics”. Thus, as in Ref. com (a), it can be used to perform “$T$-jump” simulations. These may be useful in analyzing phenomena observed in recent real $T$-jump experiments; nonequilbrium sticky bond behavior is also expected to play a role in self healing AP systems Cordier et al. (2008). Note, for example, that the timescale over which $p_{active}$ changes in Figure 11 is smaller than the equilibrium $\tau_{sb}$ ($1.5\cdot 10^{3}\tau_{LJ}$ for $h=10u_{0}$). Similarly, the timescale of self healing at a fractured surface (where $p_{active}$ is out of equilibrium) was found to be smaller than the time scale for near-equilibrium creep relaxation Cordier et al. (2008).

III.5 Nonlinear and Nonequilibrium Mechanical Properties

Figure 12a shows results for creep tests of two systems with different thermodyamics and kinetics but the same $\tau_{sb}$ ($\tau_{sb}\simeq 1.5\cdot 10^{4}\tau_{LJ}$). Both tests were performed at $k_{B}T=1.0u_{0}$. The applied stress difference $|\sigma_{z}-(\sigma_{x}+\sigma_{y})/2|=.01u_{0}/a^{3}$ is small. At times $t\ll\tau_{sb}$, the extension ratio $\lambda_{z}=L_{z}/L_{z}^{0}$ is the same for both systems. For $t\gg\tau_{sb}$, $\lambda_{z}$ is nearly linear in $ln(t)$, implying that the flow is nearly-linear creep. The system with stronger bonds and greater SB recombination shows greater resistance to flow, i. e. a smaller creep compliance.

It is interesting to relate the creep response to the quiescent dynamics. Figure 12b shows (quiescent) diffusion in the same systems. The creep response and diffusion are remarkably similar; the onset of more rapid creep in the $h=10u_{0}$, $\tau_{MC}=10\tau_{LJ}$ system under stress corresponds directly to the onset of (relative) delocalization in the quiescent state. This is consistent with a recent experiment showing connections between creep behavior and linear rheology in reversible supramolecular networks Cordier et al. (2008).

Next we consider constant volume tension simulations at $h=11.25u_{0}$, $k_{B}T=u_{0}$, and various $\tau_{MC}$. These simulations can be considered to be an extension of Ref. Rottach et al. (2007), which allowed breaking and formation of interchain bonds only at a few (discrete) strains; here SBs break and reform continuously. Here we present results for $\dot{\epsilon}=10^{-5.5}/\tau_{LJ}$. Simulations at other $\dot{\epsilon}$ were considered; larger values $\dot{\epsilon}\gtrsim t_{cage}^{-1}$ make the non-SB-related viscous stress contribution unacceptably large, while smaller values lead to more sticky bond breaking/formation during deformation than is desirable at the values of $h$ and $\tau_{MC}$ considered.

Figure 13 shows the stress difference $|\sigma_{z}-(\sigma_{x}+\sigma_{y})/2|$. With MC deactivated during deformation ($f_{MC}=0$ or equivalently $\tau_{MC}=\infty$), the stress takes a form close to that predicted by entropic elasticity Treloar (1975): $\sigma=G_{e}g(\lambda)$, where $g(\lambda)=\lambda^{2}-1/\lambda$ and $\lambda=L_{z}/L_{z}^{0}$ as above. $G_{e}$ is predicted to be $Nc_{st}p_{inter}k_{B}T/2$, where $p_{inter}<p^{*}$ is the interchain portion of active SBs; the actual value from the fit, $G_{e}=.028u_{0}/a^{3}$, is close to the predicted value $.031u_{0}/a^{3}$, indicating viscous stresses are low at this strain rate. The fit is performed for $g\leq 5$; the nonlinear behavior observed at higher $g$ arises from finite extensibility of chain segments between crosslinks (as in standard nonlinear rubber elasticity Treloar (1975); Arruda and Boyce (1993)).

The simplest result, assuming that $\dot{\epsilon}^{-1}\gg\tau_{free}$, $\tau_{sb}\gg\tau_{free}$ (here $\tau_{free}$ is the lifetime of unbonded SMs), and that SB breaking and formation rates do not vary with stress/strain, so that stress memory is lost like $exp(-t/\tau_{sb})$, is Green and Tobolsky (1946)

[TABLE]

where the first term in brackets is classical rubber elasticity and the second two terms reflect new SBs created during deformation. The $ln(\lambda)$ comes from the constant true strain rate $\dot{\epsilon}=\partial ln(\lambda)/\partial t$.

In Figure 13, stress-strain results from simulations are compared to predictions from Eq. 18 using values for $\tau_{sb}$ from Table 1 and the value of $G_{e}$ from the $\tau_{MC}=\infty$ system are shown. We confine the comparison to the linear regime ($g\leq 5$) to avoid confusion. Sticky bond recombination might be expected to slow relaxation. However, stress relaxation is actually faster than predicted by Eq. 18. We have verified that $p_{active}$ does not decrease during deformation. It appears that instead, $\tau_{sb}$ is reduced by stress. A detailed examination of this effect and comparison of nonlinear mechanical properties to theories, e. g.  Ref. ind and transient network models, e. g. Refs. Green and Tobolsky (1946); Tanaka and Edwards (1992); Vaccaro and Marruci (2000), is deferred to later work, but the data presented above suggests traditional theories will break down in the nonlinear regime.

IV Discussion and Conclusions

We have performed an initial set of simulations using a new coarse-grained model for associating polymers. The MD/MC hybrid algorithm and variable chemical kinetics allow for greater realism and flexibility than in previous simulations of AP networks. Further, the 1-1 sticky monomer binding topology imposed here reflects current experimental trends. The model was extensively validated and is able to accurately model equilibrium dynamical properties, nonlinear mechanical properties, and far-from-equilibrium systems. We studied the model over a very broad parameter space. While have emphasized that we study APs by analogy because simulations of chemically realistic AP networks are not yet computationally feasible foo (e), our results should nevertheless aid in “rational” Loveless et al. (2005) design of AP systems, especially in “transition” regimes like those discussed in this paper.

The key results presented here focused on separation, comparison and contrast of thermodynamic and chemical-kinetic effects on SB recombination, the motion of individual chains, and bulk mechanical properties. As expected, instantaneous network structure was independent of kinetics at fixed thermodynamic conditions (i. e. sticky bond strength $h/k_{B}T$), and relaxation times increases with increasing $h/k_{B}T$. Similarly, at fixed $h$ and $k_{B}T$, relaxation slows as the chemical kinetics are slowed. This was illustrated by measurements of monomer diffusion. In the physical gel regime, monomers experience a temporary “caging” similar to that found in glasses. This caging effect strengthens as SB strength is increased Kumar and Douglas (2001); Baljon et al. (2007) but also as kinetics are slowed at fixed SB strength. Analyses showed that chains become increasingly localized in a manner similar to that associated with the increase in dynamical heterogeneity in non-AP melts approaching the glass transition. Of course, the analogy should not be taken too far; in AP networks the caging is produced only by sticky monomers while in systems approaching $T_{g}$ it is produced by hard core repulsions of all monomers.

We find, as expected, that the chemical kinetics controlling $\tau_{sb}$ are “mean-field” Cates (1987) and mappable to a two-state Arrhenius model. However, as kinetic rates are increased and SB recombination becomes non-kinetically-limited, the relation between the SB lifetime $\tau_{sb}$ and other relaxation times, such as the effective SB lifetime $\tau^{*}$, becomes decidedly nontrivial. This was explicitly related to the crossover to diffusion limited SB recombination. A new quantitative relation between $\tau^{*}$ and $\tau_{sb}$ was found. Such relations should be of interest because rheological experiments can typically only access $\tau^{*}$, which is assumed to control stress relaxation (e. g., because scission followed by quick recombination does not relax stress) Cates (1987); Rotello and S. Thayumanavan (2008); Feldman et al. (2009).

While the results for $T$-variation and mechanical properties presented here were limited and somewhat preliminary, we showed examples which illustrate important effects. Analysis of diffusion on intermediate time scales illustrated the point that sticky bond and underlying polymeric timescales will in general vary differently with $T$, affecting the ‘interplay’ in nontrivial ways. In two systems with the same $\tau_{sb}$, systems with different propensities for SB recombination showed the same creep flow at short times, but those with greater recombination showed a smaller long-time (“DC”) creep compliance. These differences were directly related to the faster delocalization of chains in the quiescent state for systems with less recombination. Constant-volume deformation studies showed that, as expected, $\tau_{sb}$ is reduced by stress. Extensive studies of the variation with $T$ (in systems including attractive nonbond interactions for greater realism) and more detailed analyses of nonlinear mechanical properties are underway.

Nearly all published analytic theories for AP networks assume a single controlling relaxation time, either $\tau_{sb}$ or $\tau^{*}$, controls the ultimate relaxation properties (i. e. other relaxation times scale with the controlling time). We showed though various measurements that there is a broad parameter space (both in terms of SB strength and kinetics) within the physical gel regime where the “scaling” assumption fails. This parameter space corresponds to the conditions (1) $\tau^{*}/\tau_{sb}$ is larger but not “much larger” than unity, and/or (2) SB recombination is not kinetically limited. Deviations from this “scaling” behavior due to multiple controlling relaxation times have been observed Loveless et al. (2005); xu , but had not yet been well understood. These deviations had been previously assumed to arise from chemical disorder, and this is no doubt partially correct, but as discussed in this paper, they also arise from the ‘interplay’ between SB thermodynamics, kinetics, and polymer physics. If either (1) or (2) hold, both traditional Tanaka and Edwards (1992); Vaccaro and Marruci (2000) and more sophisticated Leibler et al. (1991); Rubinstein and Semenov (1998, 2001) theories should fail to predict the mechanical properties of AP networks. This is not meant as a criticism of the theories, merely an observation that there is a broad parameter space where one or more of their assumptions fail.

The DL and KL limits have been discussed by O’Shaughnessy and Yu O’Shaughnessy and Yu (1995); they respectively correspond to dominance of the $K$-term and $C$-term in our Eq. (10). Conditions under which systems may lie outside the KL limit and/or evidence for systems which lie outside it are also discussed, to some extent, in the context of AP networks in Refs. Loverde et al. (2005); Manassero et al. (2005); Yount et al. (2005). Coupling between SB and polymeric relaxation has also been treated approximately by Cates Cates (1987) and Leibler et. al. Leibler et al. (1991), respectively for linear EP systems and AP networks where recombination is improbable. Among published analytic theories for AP networks, Refs. Leibler et al. (1991); O’Shaughnessy and Vavylonis (1998) qualitatively treat non-kinetically-limited systems and Ref. Rubinstein and Semenov (1998) treats SB recombination. Our results are consistent with the argument of Ref. O’Shaughnessy and Vavylonis (1998) that reaction rates (here defined as non-recombinative SB exchange) reach the mean-field/KL regime only when reaction is slow compared to the longest “underlying” polymeric relaxation time [in our case $min(t_{cage},t_{chem})$] foo (k). Interestingly, Refs. de Gennes ; Friedman and O’Shaughnessy (1988) suggest that MF kinetics would apply to $\tau^{*}$ as well as $\tau_{sb}$ in dimensions $d\geq 4$ because $<|\delta\vec{r}|^{2}>\sim t^{y}$ necesarily has $y\leq 4$ (i. e. no diffusion-limited regime is possible). This suggests that diffusion-limited SB recombination, which increases $\tau^{*}/\tau_{sb}$, will increase in importance in AP systems with effectively reduced dimensionality (e. g. very thin films or “pores”). Combining the approaches of Refs. Leibler et al. (1991); O’Shaughnessy and Vavylonis (1998); Rubinstein and Semenov (1998) may be useful for developing optimal analytic theories of these systems, at least for $T$ well above $T_{g}$. However, we are not aware of any quantitative discussion of the crossover between the DL and KL regimes such as presented here foo (l).

The rheologically simple (i. e. all key relaxation times scale with $\tau_{sb}$ Yount et al. (2005); Loveless et al. (2005)) behavior observed in the majority of experiments on AP networks indicates they exhibit KL behavior. Note that these experiments have shown KL behavior even though their values of $p^{*}$ are comparable to values for systems which in our model exhibit DL behavior at low $\tau_{MC}$. This likely arises from the slow kinetics caused by the bulkiness and directional interactions of real sticky monomers. Creating (real) strong-binding SMs with even faster kinetics seems to be difficult. However, one can move out of the KL regime (at fixed kinetic rates) simply by slowing the polymeric relaxation times, e. g. by going to higher concentrations and/or entangled chains. Experiments in this regime are underway xu , and seem to show a breakdown of the simple scaling; for example, they show an unusually high power law dependence of viscosity on concentration, which appears to arise because $\tau_{sb}$ is of order the time scale for reptation. Ref. Loveless et al. (2005) also shows an apparent (if weak) breakdown in scaling at the highest frequencies considered. In the context of these observations, we note that the parameter space where (1) and or (2) hold may be of greatest interest for designing materials with novel mechanical properties. Our model seems well suited to aid in understanding the complicated behavior of AP networks in this regime.

Here we have left the regime of physically entangled APs, which is the regime treated by some key analytic AP theories Leibler et al. (1991); Rubinstein and Semenov (2001); Semenov (2006), untouched. Also, the “interplay” described in this paper should depend on the details of sticky monomer arrangement along chains, not just $N$ and $c_{st}$. Studies of systems with a wide range of $N$ and $c_{st}$, as well as inhomogeneous (chemically disordered) systems, are underway.

In real AP networks the sticky and regular monomers have different sizes and chemistries. Thus an obvious extension of our model would be to increase the differences between the sticky and normal monomers. For example, changing secondary interactions may induce microphase separation Manassero et al. (2005); Grest et al. (1996). Another extension would be use of a more realistic sticky bonding potential such as those used to model H-bonds (see e. g.  Ref. kar and refs. therein), but here we have focused on chemistry-independent properties. A coarse-grained way to to capture this would be to keep the binary bonding rules, but include more than one sticky site per SM; this would increase the directionality of bonding, which is a key to the performance of the UPy systems Sijbesma et al. (1997). Finally, nanocomposites of associating polymers Sprakel et al. (2007), where the presence of nanoparticles may or may not Sen et al. (2007); Allegra et al. (2008) affect single-chain structure but will certainly affect AP network structure, should have even richer physics than regular polymer nanocomposites.

Michael Rubinstein, Kathleen E. Feldman, Arlette R. C. Baljon, Phillip A. Pincus, Edward J. Kramer, Frank L. H. Brown, Stephen L. Craig, Gary S. Grest, Richard C. Elliott, Stuart J. Rowan, Steven J. Plimpton, and Ronald G. Larson provided helpful discussions. Michael Rubinstein additionally provided helpful comments on the manuscript. Brian B. Rochford assisted with the numerical analysis of $P_{recomb}$. This work was supported by the MRSEC Program of the National Science Foundation under Award No. DMR05-20415.

Appendix A Percolation Gel Transition

Let $N^{cl}_{i}\equiv N_{ch}P(N_{ch},i)$ be the average number of disconnected clusters of $i$ chains in a system of $N_{ch}$ chains with periodic boundary conditions (at any given time). $P(N_{ch},i)$ is a cluster size probability distribution with $\sum_{i=1}^{N_{ch}}iP(N_{ch},i)\equiv 1$. The number averaged cluster size is then $N_{n}=N_{ch}\sum_{i=1}^{N_{ch}}iP(N_{ch},i)$ and the weight averaged cluster size is $N_{w}=N_{n}^{-1}N_{ch}\sum_{i=1}^{N_{ch}}i^{2}P(N_{ch},i)$. For an infinite system, the percolation gel transition occurs (by definition) when $p^{*}$ exceeds $p_{perc}$; $N_{w}$ diverges at $p^{*}=p_{perc}$ Rubinstein and Colby (2003). However, computer simulations are limited to finite $N_{ch}$, and the value of $N_{w}$ cannot exceed $N_{ch}$. Thus the $N_{ch}$-dependent geometric percolation $p^{*}=p_{span}$ (at which one aggegrate spans the system) approaches $p_{perc}$ from below as $N_{ch}\to\infty$) com (c). Fortunately, $p_{perc}$ and $h_{perc}$ (the value of $h$ at which $p^{*}=p_{perc}$ in an infinite system) can be estimated for finite $N_{ch}$ using a standard finite size analysis Frenkel and Smit (2002).

We perform such an analysis, following Ref. Pandey and Landau (1985). Figure 14 shows this analysis for $N=50$, $c_{st}=.08$ systems at $k_{B}T=1.0u_{0}$. The figure plots the rescaled variables $(N_{w}/N_{ch})^{-\gamma}N_{ch}^{-\gamma/3\nu}$ vs. $((p_{perc}-p^{*})/p^{*})N_{ch}^{1/3\nu}$ Pandey and Landau (1985). $p^{perc}\simeq.40$ is close to the predicted value $1/(Nc_{st}p^{*}f_{inter}-1)$ com (c), where $f_{inter}$ is the fraction of SBs which are interchain rather than intrachain. The exponents used to collapse the data in the figure are $\gamma\simeq 1.7$ and $\nu\simeq 1.2$; considering a narrower range of $h$ gives values consistent with predictions from the theory of critical phenomena ($\gamma\simeq 1.8,\ \nu\simeq 0.9$) de Gennes (1979). These exponents been extensively discussed in the literature Pandey and Landau (1985); com (c); Y.Liu and Pandey (1996) and need not be discussed further here. In the figure, $h$ increases going from right to left, and percolation occurs at $h_{perc}\simeq 4.25u_{0}$.

Appendix B Numerical Analysis of SB Recombination

In this paper, the sticky bond self-correlation function $P_{auto}(\triangle t)$ is the probability that a bond between two given SMs exists both at times $t$ and $t+\triangle t$, while the SB “transition function” $P_{trans}(\triangle t)$ is the probability that the bond exists continuously between times $t$ and $t+\triangle t$. The SB “recombination function” $P_{recomb}(\triangle t)\equiv P_{auto}(\triangle t)-P_{trans}(\triangle t)$ is the probability that a pair of SMs will be bonded at two times separated by $\triangle t$ but that the bond between them has broken at least once during that interval. We find $P_{trans}$ exhibits nearly single-exponential decay, $P_{trans}=\exp(-t/\tau_{sb})$ for all systems, while for $h\gtrsim 10u_{0}$, $P_{auto}$ also shows exponential decay for (at least) the first decade. Values of $\tau^{*}$ presented here are measured from fits to $P_{trans}=\exp(-t/\tau^{*})$. All quantities are averaged over all SM pairs and all $t$.