Irreversible Step-Growth Polymerization with weighted dynamic networks

Many research fields benefit from recent developments in random graphs and network theory. In particular, polymerization process can be described well with random graph models [8, 6, 3] and these results are consistent with well-known models, e.g. Flory [1] or Stockmayer [10].

Conventional polymer networks are formed by a process called polymerization during which monomers bind together by means of covalent bonds and form large molecular structures. One of the most common polymerization processes is the step-growth polymerization of multifunctional monomers. This process leads to hyperbranched polymers of irregular topologies that undergo a phase transition in their connectivity structure which is marked by emergence of the gel, i.e. the giant molecule that spans the whole volume. Flory provided simple analytical expressions for the average molecular size for several particular configurations. Later Stockmayer presented a formal expression for the whole distribution of molecular sizes, but its application is quite limited because of complicated computations.

Here we should mention fast approximated methods as in papers [12, 11, 2, 7]. However, these are very specific and hard to adapt to new polymerization schemes. Molecular dynamic simulation approaches provide lots of details on the structure of polymer networks, however they are computationally expensive and, thus, limited to small samples and short time scales.

On the other hand random graph approaches present a opportunity to obtain exactly solvable expressions for hyperbranched polymers. In the case of symmetric covalent bonds the structure of polymer networks is described by the configuration model for undirected random networks [8, 6]. Developments in directed configuration models allowed to develop a generic polymerization framework that covers asymmetrical bonds as well.

Unfortunately, these models miss the case when both symmetrical and asymmetrical bonds are present in the system, and do not account for rate at which these bonds are formed. This paper tries to address these issues.

We start from considering simple case when we have two types of monomers, say, A and B. A can connect with B, but not with self; on the opposite, B can connect both with A and B. This differs from the models considered in other papers by mixing both directional and non-directional connections.

We distinguish monomer species by counting the numbers of functional groups of both types $I$, $J$ and the numbers of in- and bidirectional-edges $i,j$. During the progress of polymerisation the functional groups are converted into chemical bonds between the monomers, and the concentration profiles $M_{i,j,I,J}(t)$ obey the following master equation:

$$\begin{array}[]{l}\frac{\partial}{\partial t}M_{i,j,I,J}(t)=(I-i+1)(\nu_{01}-\mu_{01}(t))M_{i-1,j,I,J}(t)+\\ (J-j+1)(\nu_{10}-\mu_{10}(t)+\nu_{01}-\mu_{01}(t))M_{i,j-1,I,J}(t)-\\ ((I-i)(\nu_{01}-\mu_{01}(t))+(J-j)(\nu_{10}-\mu_{10}(t)+\nu_{01}-\mu_{01}(t)))M_{i,j,I,J}(t)\end{array}$$

(1)

Like in the case of fully oriented graph, considered in [3], probability of connecting functional group A is proportional to the number of free functional groups B. On the other hand, since the functional groups B can connect both to A and B, probability of connecting them is proportional to the total number of free functional groups.

Initially at $t=0$ there are no bonds, thus we have the following initial condition

$$\begin{array}[]{ll}M_{i,j,I,J}(0)=P(I,J)&\mbox{if }i=j=0,\\ M_{i,j,I,J}(0)=0&\mbox{otherwise.}\end{array}$$

(2)

Here $P(I,J)$ denotes initial distribution of monomers of $I,J$ type. Define $\nu_{mn}$ to be partial moments of the initial monomer type distribution

$$\nu_{mn}=\sum_{I,J\geq 0}I^{m}J^{n}P(I,J)$$

Let the time dependent degree distribution $u(i,j,t)=\sum_{I,J}M_{i,j,I,J}(t)$, then

$$\mu_{mn}(t)=\sum_{i,j\geq 0}i^{m}j^{n}u(i,j,t),$$

(3)

e.g. $\mu_{mn}$ are partial moments of degree distribution. Note that above definition is not strict since summation is performed up to $I$ and $J$ correspondingly, so one should substitute expression for $u(i,j,t)$ and change summation order, e.g. $\mu_{mn}(t)=\sum_{I,J}\sum_{i,j}i^{m}j^{n}M_{i,j,I,J}(t)$.

In order to solve the system of differential equations (1,2) we proceed to generating functions, namely, let $m(x,y,I,J,t)=\sum_{i=0}^{I}\sum_{j=0}^{J}x^{i}y^{j}M_{i,j,I,J}(t)$.

Obviously, we have the following relations:

$$\begin{array}[]{l}\frac{\partial}{\partial x}m(x,y,I,J,t)=\sum_{i\geq 1,j}ix^{i-1}y^{j}M_{i,j,I,J}(t),\\ \frac{\partial}{\partial y}m(x,y,I,J,t)=\sum_{j\geq 1,i}jx^{i}y^{j-1}M_{i,j,I,J}(t),\\ m(x,y,I,J,0)=P(I,J).\end{array}$$

(4)

Let us multiply equations (1) by $x^{i}y^{j}$ and sum up over all $i,j$. Taking advantage of (4), we obtain the following system of PDEs

$$\begin{array}[]{l}\frac{\partial m(x,y,I,J,t)}{\partial t}=(\nu_{01}-\mu_{01}(t))(Ixm(x,y,I,J,t)-x^{2}\frac{\partial m(x,y,I,J,t)}{\partial x})+\\ (\nu_{10}-\mu_{10}(t)+\nu_{01}-\mu_{01}(t))(Jym(x,y,I,J,t)-y^{2}\frac{\partial m(x,y,I,J,t)}{\partial y})-\\ (\nu_{01}-\mu_{01}(t))(Im(x,y,I,J,t)-x\frac{\partial m(x,y,I,J,t)}{\partial x})-\\ (\nu_{01}-\mu_{01}(t)+\mu_{10}-\mu_{10}(t))(J(x,y,I,J,t)m-y\frac{\partial m(x,y,I,J,t)}{y}),\\ m(x,y,I,J,0)=P(I,J).\end{array}$$

(5)

Since $u(i,j,t)=\sum_{I,J}M_{i,j,I,J}(t)$, we can write

$$\begin{array}[]{l}\mu_{10}(t)=\sum_{I,J}\left.\frac{\partial m(x,y,I,J,t)}{\partial x}\right|_{x=1,y=1},\\ \mu_{01}(t)=\sum_{I,J}\left.\frac{\partial m(x,y,I,J,t)}{\partial y}\right|_{x=1,y=1}.\end{array}$$

(6)

Next, we’re going to obtain a system of ODEs for the functions $\mu_{10}(t)$ and $\mu_{01}(t)$.

Let us rewrite the differential equation from (5) as

$$\begin{array}[]{l}\frac{\partial m(x,y,I,J,t)}{\partial t}=(\nu_{01}-\mu_{01}(t))(x-1)(Im(x,y,I,J,t)-x\frac{\partial m(x,y,I,J,t)}{\partial x})+\\ (\nu_{01}-\mu_{10}(t)+\mu_{01}-\mu_{01}(t))(y-1)(Jm(x,y,I,J,t)-y\frac{\partial m(x,y,I,J,t)}{\partial y})\end{array}$$

(7)

Differentiating (7) by $x$ and taking $x=1,y=1$, we have

$$\left.\frac{\partial m(x,y,I,J,t)}{\partial t\partial x}\right|_{x=y=1}=(\nu_{01}-\mu_{01}(t))\left(Im(1,1,I,J,t)-\left.\frac{\partial m(x,y,I,J,t)}{\partial x}\right|_{x=y=1}\right)$$

After summing up over all $I,J$ and taking into account (6) we have

$$\frac{\partial}{\partial t}\mu_{10}(t)=(\nu_{01}-\mu_{01}(t))(\sum_{I,J}Im(1,1,I,J,t)-\mu_{10}(t)).$$

Note that $\sum_{I,J}Im(1,1,I,J,t)=\sum_{I,J}I\sum_{i,j}M_{i,j,I,J}(t)$. Since total number of monomer species of each type remains the same, $\sum_{i,j}M_{i,j,I,J}(t)=\sum_{i,j}M_{i,j,I,J}(0)$, thus $\sum_{I,J}Im(1,1,I,J,t)=\sum_{I,J}IP(I,J)=\nu_{10}$.

Finally, $\mu_{10}(0)=\sum_{I,J}\left.\frac{\partial m(x,y,I,J,t)}{\partial x}\right|_{x=y=1,t=0}$ and it remains to recall that $m(x,y,I,J,0)=P(I,J)$ for all $x,y$, so $\mu_{10}(0)=0$.

Differentiating (7) by $y$, taking $x=1,y=1$, summing up by $I,J$ and taking into account (6), we obtain system of ODEs for $\mu_{01}(t)$ and $\mu_{10}(t)$.

$$\begin{array}[]{l}\frac{\partial\mu_{10}(t)}{\partial t}=(\nu_{01}-\mu_{01}(t))(\nu_{10}-\mu_{10}(t)),\\ \frac{\partial\mu_{01}(t)}{\partial t}=(\nu_{01}-\mu_{01}(t)+\nu_{10}-\mu_{10}(t))(\nu_{01}-\mu_{01}(t)),\\ \mu_{01}(0)=\mu_{10}(0)=0.\end{array}$$

(8)

In [4] they consider fully oriented graph, which makes the system (8) symmetrical, e.g. the second equation looks like $\frac{\partial\mu_{01}(t)}{t}=(\nu_{10}-\mu_{10}(t))(\nu_{01}-\mu_{01}(t))$; evidently, one can conclude $\mu_{01}(t)=\mu_{10}(t)$ and can easily find analytical solution of such a system.

In our case $\mu_{01}(t)\not\equiv\mu_{10}(t)$, so we’re solving (8) numerically with RK methods.

Now, let’s go back to solving the system (5). Assume solution in form

$$m(x,y,I,J,t)=\left(\frac{1+A(t)t}{1+A(t)tx}\right)^{-I}\left(\frac{1+B(t)t}{1+B(t)ty}\right)^{-J}P(I,J).$$

Consequently,

$$\begin{array}[]{l}\log m(x,y,I,J,t)=-I\log(1+A(t)t)+I\log(1+A(t)tx)-\\ J\log(1+B(t)t)+J\log(1+B(t)ty)+\log P(I,J).\end{array}$$

(9)

Let us rewrite (7), dividing left and right parts by $m(x,y,I,J,t)$:

$$\begin{array}[]{l}\frac{\partial\log m(x,y,I,J,t)}{\partial t}=(\nu_{01}-\mu_{01}(t))(x-1)(I-x\frac{\partial\log m(x,y,I,J,t)}{\partial x})+\\ (\nu_{10}-\mu_{10}(t)+\nu_{01}-\mu_{01}(t))(y-1)(J-y\frac{\partial\log m(x,y,I,J,t)}{\partial y})\end{array}$$

and substitute (9).

Equate the terms depending on $A(t)$:

$$\frac{I(A(t)+t\frac{\partial A(t)}{\partial t})(x-1)}{(1+A(t)t)(1+A(t)tx)}=\frac{I(\nu_{01}-\mu_{01}(t))(x-1)}{1+A(t)tx}$$

Thus,

$$\frac{A(t)+t\frac{\partial A(t)}{\partial t}}{1+A(t)t}=\nu_{01}-\mu_{01}(t)$$

with initial condition is $A(0)=\nu_{01}$. Consequently,

$$A(t)=\frac{\exp\{\int_{0}^{t}\nu_{01}-\mu_{01}(u)du\}-1}{t}.$$

(10)

Similarly,

$$B(t)=\frac{\exp\{\int_{0}^{t}\nu_{10}-\mu_{10}(u)+\nu_{01}-\mu_{01}(u)du\}-1}{t}.$$

(11)

Since $\sum_{i=0}^{n}\binom{n}{k}x^{i}p^{i}(1+p)^{-n}=\left(\frac{1+p}{1+xp}\right)^{-n}$, we get expression for $M_{i,j,I,J}(t)$:

$$M_{i,j,I,J}(t)=\binom{I}{i}(tA(t))^{i}(1+tA(t))^{-I}\binom{J}{j}(tB(t))^{j}(1+tB(t))^{-J}P(I,J).$$

This expression can be rewritten using binominal probabilities. Let us denote $p_{A}(t)=\frac{tA(t)}{1+tA(t)}$ and $p_{B}(t)=\frac{tB(t)}{1+tB(t)}$, then

$$M_{i,j,I,J}(t)=\binom{I}{i}p_{A}(t)^{i}(1-p_{A}(t))^{I-i}\binom{J}{j}p_{B}(t)^{j}(1-p_{B}(t))^{J-j}P(I,J).$$

(12)

Consequently,

$$u(i,j,t)=\sum_{I,J}\binom{I}{i}p_{A}(t)^{i}(1-p_{A}(t))^{I-i}\binom{J}{j}p_{B}(t)^{j}(1-p_{B}(t))^{J-j}P(I,J)$$

(13)

Once expression for $u(i,j,t)$ is set up, we can calculate $\mu_{nm}(t)$ according to (3), namely

$$\begin{array}[]{l}\mu_{00}(t)=1,\\ \mu_{10}(t)=p_{A}(t)\nu_{10},\mu_{01}(t)=p_{B}(t)\nu_{01},\\ \mu_{20}(t)=p_{A}(t)(p_{A}(t)\nu_{20}-p_{A}(t)\nu_{10}+\nu_{10}),\\ \mu_{02}(t)=p_{B}(t)(p_{B}(t)\nu_{02}-p_{B}(t)\nu_{01}+\nu_{01}),\\ \mu_{11}(t)=p_{A}(t)p_{B}(t)\nu_{11},\mu_{22}(t)=p_{A}(t)p_{B}(t)(p_{A}(t)p_{B}(t)(\nu_{22}-\nu_{21}-\nu_{12}+\nu_{11})+\\ p_{A}(t)(\nu_{21}-\nu_{11})+p_{B}(t)(\nu_{12}-\nu_{11})+\nu_{11}).\end{array}$$

(14)

Note that obtained relations for $\mu_{10}(t)$ and $\mu_{01}(t)$ allow much faster calculation of $p_{A}(t)$ and $p_{B}(t)$ compared to their definitions. Stationary solution for the system (8) can be partially obtained by equating left-hand sides to zero, thus we immediately have $\mu_{01}(t)\rightarrow\nu_{01}$ as $t\rightarrow\infty$.

Systems with a single monomer type

Consider the case when there is only one monomer species bearing $I$ groups of type $A$ and $J$ groups of type $B$. This corresponds to condition $P(I,J)=1$, thus all the sums over all $I$ and $J$ reduce to a single summand. Moments of distribution $P$ are $\nu_{mn}=I^{m}J^{n}$.

Unfortunatelly, even in this case the system (8) cannot be solved analytically, and we should stick to numerical solution.

Consequently we can calculate $A(t),B(t)$, and $p_{A}(t),p_{B}(t)$. Numerical results for this kind of model are presented in Appendix.

Next, we’ll try to derive gen transition point. Let $U(x,y,t)$ be generating function for degree distribution $u(i,j,t)$. Suppose, we select a vertex that is at the end of a directed randomly-chosen edge. Conditional degree distribution in this case is $u_{in}(n,k,t)=\frac{n}{\mu_{10}(t)}u(n,k,t)$. Correspondingly, if we select a vertex at non-directed randomly-chosen edge, conditional degree distribution would be $u_{bi}(n,k,t)=\frac{k}{\mu_{01}(t)}u(n,k,t)$. Corresponding generating functions can be expressed as

$$\begin{array}[]{l}U_{in}(x,y,t)=\frac{1}{\mu_{10}(t)}\frac{\partial U(x,y,t)}{\partial x},\\ U_{bi}(x,y,t)=\frac{1}{\mu_{01}(t)}\frac{\partial U(x,y,t)}{\partial y}.\end{array}$$

(15)

If $w(s,t)$ is the probability that a randomly chosen monomer belongs to a component of size $s$ (e.g. molecular weight distribution), then $w(1,t)=u(0,0,t),w(2,t)=u(1,0,t)u(0,1,t)\left(\frac{1}{\mu_{01}(t)}+\frac{1}{\mu_{10}(t)}\right)+\frac{u(0,1,t)^{2}}{\mu_{01}(t)}$, etc.

Let $W(x,t)$ be GF of $w(s,t)$. Now consider biased choice for the starting vertex. Suppose, one selects a directed edge at random and the end vertex of this edge as root. Then, $w_{in}(s,t)$ denotes the distribution of weak-component sizes associated with this root; corresponding GF we denote through $W_{in}(x,t)$. Similarly, for the non-directed edge we define GF $W_{bi}(x,t)$. Bind $W_{in},W_{bi},W,U_{in},U_{bi},U$ together.

Let us start by selecting a root vertex that we arrive at by following a directed edge. The probability of the root to have $n$ in-edges and $k$ non-directed edges is $u_{in}(n,k,t)$. Each of the in-edges is associated with a weak component of the size $w_{bi}(s,t)$, thus the sum of sizes of all components reached through the in-edges is distributed as a sum of independent copies of $w_{bi}(s,t)$. On the other hand, each non-directional edge can be associated with a weak component of the size $w_{in}(s,t)$ with conditional probability proportional to $\mu_{10}(t)$ and with a weak component of the size $w_{bi}(s,t)$ with probability proportional to $\mu_{01}(t)$. Let $\kappa(t)=\frac{\mu_{10}(t)}{\mu_{10}(t)+\mu_{01}(t)}$. Introduce $W_{+}(s,t)=\kappa(t)W_{in}(x,t)+(1-\kappa(t))W_{bi}(x,t)$, then the generating function for a weak component distrubution reached through the non-directional edge is $W_{+}(s,t)$.

Thus, the distribution for the sum of sizes of all components originated at the root is obtained as

$$\sum_{n,k}u_{in}(n,k,t)W_{bi}(x,t)^{n}W_{+}(x,t)^{k}=U_{in}(W_{bi}(x,t),W_{+}(x,t)).$$

On the other hand, the total number of all vertices reachable from the root plus one can be also considered as the size of the weak component reached by following the original in-edge. Thus, one obtains the following relation:

$$W_{in}(x)=xU_{in}(W_{bi}(x,t),W_{+}(x,t)).$$

Similarly, we get the following system of functional equations:

$$\begin{array}[]{l}W(x,t)=xU(W_{bi}(x,t),W_{+}(x,t),t),\\ W_{in}(x,t)=xU_{in}(W_{bi}(x,t),W_{+}(x,t),t),\\ W_{bi}(x,t)=xU_{bi}(W_{bi}(x,t),W_{+}(x,t),t),\end{array}$$

or equivalently,

$$\begin{array}[]{l}W(x,t)=xU(W_{bi}(x,t),W_{+}(x,t),t),\\ W_{+}(x,t)=xU_{+}(W_{bi}(x,t),W_{+}(x,t),t),\\ W_{bi}(x,t)=xU_{bi}(W_{bi}(x,t),W_{+}(x,t),t),\end{array}$$

(16)

where $U_{+}(n,k,t)=\kappa(t)U_{in}(n,k,t)+(1-\kappa)U_{bi}(n,k,t)$.

Weight-average molecular weight is calculated as $\frac{W^{\prime}(1,t)}{W(1)}$.

Since $\mu_{mn}(t)=\left.\left(x\frac{\partial}{\partial x}\right)^{m}\left(y\frac{\partial}{\partial y}\right)^{n}U(x,y,t)\right|_{x=y=1}$, we have $\mu_{10}(t)=U_{in}(1,1,t)\mu_{10}(t)$, e.g. $U_{in}(1,1,t)=1$. Similarly we have $U_{bi}(1,1,t)=1$.

Consequently, (16) gives us $W_{in}(1,t)=1$ and $W_{bi}(1,t)=1$. It remains to calculate derivative of $W(x,t)$. Let us take derivative in $x$ for all relations in (16) at $x=1$:

$$\begin{array}[]{l}\frac{\partial W(1,t)}{\partial x}=1+\frac{\partial U(1,1,t)}{\partial x}\frac{\partial W_{bi}(1,t)}{\partial x}+\frac{\partial U(1,1,t)}{\partial y}\frac{\partial W_{+}(1,t)}{\partial x}\\ \frac{\partial W_{+}(1,t)}{\partial x}=1+\frac{\partial U_{+}(1,1,t)}{\partial x}\frac{\partial W_{bi}(1,t)}{\partial x}+\frac{\partial U_{+}(1,1,t)}{\partial y}\frac{\partial W_{+}(1,t)}{\partial x}\\ \frac{\partial W_{bi}(1,t)}{\partial x}=1+\frac{\partial U_{bi}(1,1,t)}{\partial x}\frac{\partial W_{bi}(1,t)}{\partial x}+\frac{\partial U_{bi}(1,1,t)}{\partial y}\frac{\partial W_{+}(1,t)}{\partial x}.\end{array}$$

Solution to this linear system degenerates if

$$\frac{\partial U_{+}(1,1,t)}{\partial x}\frac{\partial U_{bi}(1,1,t)}{\partial y}=(\frac{\partial U_{+}(1,1,t)}{\partial y}-1)(\frac{\partial U_{bi}(1,1,t)}{\partial x}-1)$$

(17)

and corresponds to unlimited growth of the weak component.

From (15) we obtain

$$\begin{array}[]{l}\frac{\partial U_{bi}(1,1,t)}{\partial y}=\frac{\mu_{02}(t)}{\mu_{01}(t)}-1,\\ \frac{\partial U_{bi}(1,1,t)}{\partial x}=\frac{\mu_{11}(t)}{\mu_{01}(t)},\\ \frac{\partial U_{in}(1,1,t)}{\partial x}=\frac{\mu_{20}(t)}{\mu_{10}(t)}-1,\\ \frac{\partial U_{in}(1,1,t)}{\partial y}=\frac{\mu_{11}(t)}{\mu_{10}(t)}.\end{array}$$

(18)

Substituting this in its turn into the relation (17), we get the criteria for gel transition point.

We’ll use this statement to obtain numerical results in the Appendix.

Let us generalize the considered model by introducing weights which express reaction rate. Namely, let us have two types of monomers, $A$ and $B$. As before, $A$ can connect with $B$, but not with self; however, $B$ can react not only with $A$, but also with $B$, and let us account that the latter is less probable. We will express this change by introducing the weight $0\leq w\leq 1$ which denotes “likelyness” of establishing connection $B-B$.

First of all, let us see how the original master equation (1) changes:

$$\begin{array}[]{l}\frac{\partial}{\partial t}M_{i,j,I,J}(t)=(I-i+1)(\nu_{01}-\mu_{01}(t))M_{i-1,j,I,J}(t)+\\ (J-j+1)(\nu_{10}-\mu_{10}(t)+w(\nu_{01}-\mu_{01}(t)))M_{i,j-1,I,J}(t)-\\ ((I-i)(\nu_{01}-\mu_{01}(t))+(J-j)(\nu_{10}-\mu_{10}(t)+w(\nu_{01}-\mu_{01}(t))))M_{i,j,I,J}(t)\end{array}$$

(20)

while initial conditions (2) remain the same. Consequently, we proceed to GFs and obtain the following system of ODEs which resembles (8):

$$\begin{array}[]{l}\frac{\partial\mu_{10}(t)}{t}=(\nu_{01}-\mu_{01}(t))(\nu_{10}-\mu_{10}(t)),\\ \frac{\partial\mu_{01}(t)}{t}=(w(\nu_{01}-\mu_{01}(t))+\nu_{10}-\mu_{10}(t))(\nu_{01}-\mu_{01}(t)),\\ \mu_{01}(0)=\mu_{10}(0)=0.\end{array}$$

(21)