Polymerization of Branched Polyetherimides: Comparison between Monte Carlo Simulation and Flory-Stockmayer Theory

TL;DR

This study develops a Monte Carlo simulation method to model the polymerization of branched polyetherimides and compares its predictions with Flory-Stockmayer theory, highlighting the conditions under which each method is accurate.

Contribution

The paper introduces a Monte Carlo simulation approach based on the Gillespie algorithm for branched polyetherimide polymerization and compares its accuracy to Flory-Stockmayer theory across different system states.

Findings

MC simulation accurately predicts molecular weight distribution below gel point.

Flory-Stockmayer theory fails near or above gel point.

MC simulations are computationally efficient for small systems with large-scale implications.

Abstract

A Monte Carlo (MC) simulation method based on the Gillespie algorithm is developed for the polymerization of branched polyetherimides from two back-bone monomers [4,4'-bisphenol A dianhydride (BPADA) and m-phenylenediamine (MPD)], a chain terminator [phthalic anhydride (PA)], and a branching agent [tris[4-(4-aminophenoxy)phenyl] ethane (TAPE)]. This polymerization involves 4 reactions that all can be reduced to a condensation reaction between an amine group and a carboxylic anhydride group. By comparing the MC simulation results to the predictions of the Flory-Stockmayer theory on the molecular weight distribution, we show that the rates of the 4 reactions in the MC simulations should be set based on the concentrations of the functional groups on the monomers involved in each reaction. Using the MC simulations, we show that the Flory-Stockmayer theory predicts the molecular weight distribution well for systems below the gel point that is set by the functionality of the branching agent but starts to fail for systems around or above the gel point. However, for all the systems the MC simulations can be used to reliably predict the molecular weight distribution. Even for a macroscopic system, a converging distribution can be quickly obtained through the MC simulations on a system of only a few hundred to a few thousand monomers but with the same molar ratios of monomers as in the macroscopic system.

Figures

Tables

Table 1. Table 1: System used for checking the way to compute the reaction rates.

Monomer	PA	BPADA	MPD	TAPE
Number	2000	0	0	1000

Table 2. Table 2: Three fully reacted, stoichiometric systems below, around, and beyond the gel point. The first column is the system label. The next 4 columns list the number of each type of monomers. The values of ρ 1 subscript 𝜌 1 \rho_{1} and ρ 2 subscript 𝜌 2 \rho_{2} are determined from the monomer numbers. The value of α 𝛼 \alpha is computed using Eq. ( 20 ). The average molecular weights, M n subscript 𝑀 𝑛 M_{n} , M w subscript 𝑀 𝑤 M_{w} , and M z subscript 𝑀 𝑧 M_{z} , are from the MC simulations.

	PA	BPADA	MPD	TAPE	${\rho }_1$	${\rho }_2$	$p_A$	$p_B$	$\alpha$	$M_n$ (Da)	$M_w$ (Da)	$M_z$ (Da)
	50	670	680	10	0.0360	0.0216	1	1	0.366	$19120\pm 33$	$52126\pm 545$	$78671\pm 2488$
$S_{\simeq }$	50	670	671	16	0.0360	0.0345	1	1	0.481	$22000\pm 15$	$77353\pm 307$	$116227\pm 499$
$S_{>}$	50	670	620	50	0.0360	0.108	1	1	0.743	$51055\pm 125$	$330124\pm 458$	$369588\pm 348$

Table 3. Table 3: Four fully reacted, stoichiometric systems all below the gel point but with the size increased proportionally without changing the fraction of each type of monomers. The entries have the same format as in Table  2 . The subscript of the system label in the first column indicates the size ratio with respect to the base system, S 1 subscript 𝑆 1 S_{1} . The average molecular weights, M n subscript 𝑀 𝑛 M_{n} , M w subscript 𝑀 𝑤 M_{w} , and M z subscript 𝑀 𝑧 M_{z} , are from the MC simulations.

	PA	BPADA	MPD	TAPE	${\rho }_1$	${\rho }_2$	$p_A$	$p_B$	$\alpha$	$M_n$ (Da)	$M_w$ (Da)	$M_z$ (Da)
$S_1$	10	134	136	2	0.0360	0.0216	1	1	0.366	$15742\pm 14$	$30829\pm 46$	$37334\pm 57$
$S_{10}$	100	1340	1360	20	0.0360	0.0216	1	1	0.366	$19799\pm 18$	$59940\pm 607$	$101321\pm 1441$
$S_{50}$	500	6700	6800	100	0.0360	0.0216	1	1	0.366	$20361\pm 4$	$73582\pm 619$	$161904\pm 2612$
$S_{80}$	800	10720	10880	160	0.0360	0.0216	1	1	0.366	$20417\pm 3$	$75980\pm 550$	$177919\pm 2712$

Table 4. Table 4: Five partially reacted, stiochiometric systems (i.e., p A = p B < 1 subscript 𝑝 𝐴 subscript 𝑝 𝐵 1 p_{A}=p_{B}<1 ). The entries have the same format as in Table  2 . The superscript of the system label indicates the values of p A subscript 𝑝 𝐴 p_{A} and p B subscript 𝑝 𝐵 p_{B} . The first two are below and the rest three are beyond the gel point. The average molecular weights, M n subscript 𝑀 𝑛 M_{n} , M w subscript 𝑀 𝑤 M_{w} , and M z subscript 𝑀 𝑧 M_{z} , are from the MC simulations.

	PA	BPADA	MPD	TAPE	${\rho }_1$	${\rho }_2$	$p_A$	$p_B$	$\alpha$	$M_n$ (Da)	$M_w$ (Da)	$M_z$ (Da)
$S^{0.95}$	50	670	620	50	0.0360	0.108	0.95	0.95	0.419	$5965\pm 3$	$24779\pm 331$	$47952\pm 824$
$S^{0.96}$	50	670	620	50	0.0360	0.108	0.96	0.96	0.462	$7386\pm 5$	$36488\pm 488$	$69368\pm 1097$
$S^{0.97}$	50	670	620	50	0.0360	0.108	0.97	0.97	0.513	$9655\pm 10$	$58237\pm 807$	$103512\pm 1504$
$S^{0.98}$	50	670	620	50	0.0360	0.108	0.98	0.98	0.574	$13780\pm 23$	$105576\pm 1435$	$165690\pm 2080$
$S^{0.99}$	50	670	620	50	0.0360	0.108	0.99	0.99	0.649	$22710\pm 70$	$200514\pm 1915$	$266624\pm 2083$

Table 5. Table 5: Three partially reacted, nonstiochiometric systems (i.e., p A ≠ p B subscript 𝑝 𝐴 subscript 𝑝 𝐵 p_{A}\neq p_{B} and both are less than 1) below, around, and beyond the gel point. The entries have the same format as in Table  2 . The average molecular weights, M n subscript 𝑀 𝑛 M_{n} , M w subscript 𝑀 𝑤 M_{w} , and M z subscript 𝑀 𝑧 M_{z} , are from the MC simulations.

	PA	BPADA	MPD	TAPE	${\rho }_1$	${\rho }_2$	$p_A$	$p_B$	$\alpha$	$M_n$ (Da)	$M_w$ (Da)	$M_z$ (Da)
	50	670	664	50	0.0360	0.101	0.99	0.93	0.445	$7225\pm 2$	$33130\pm 195$	$62388\pm 448$
$S_{\simeq }^n$	50	670	649	50	0.0360	0.104	0.99	0.95	0.502	$9530\pm 4$	$54269\pm 349$	$96424\pm 656$
$S_{>}^n$	50	670	634	50	0.0360	0.106	0.99	0.97	0.569	$13824\pm 10$	$103074\pm 624$	$162837\pm 915$

Equations

$$p_A\sum _{q=1}^i{f}_q{A}_q=p_B\sum _{h=1}^j{g}_h{B}_h.$$

$$\begin{array}{rl}N\{m,n\}& =K\frac{\left({\sum }_{q=1}^i{f}_q{m}_q-{\sum }_{q=1}^i{m}_q\right)!}{\left({\sum }_{q=1}^i{f}_q{m}_q-{\sum }_{q=1}^i{m}_q-{\sum }_{h=1}^j{n}_h+1\right)!}\\ & \times \frac{\left({\sum }_{h=1}^j{g}_h{n}_h-{\sum }_{h=1}^j{n}_h\right)!}{\left({\sum }_{h=1}^j{g}_h{n}_h-{\sum }_{h=1}^j{n}_h-{\sum }_{q=1}^i{m}_q+1\right)!}\\ & \times \prod _{q=1}^i\frac{x_q^{m_q}}{m_q!}\prod _{h=1}^j\frac{y_h^{n_h}}{n_h!}\end{array}$$

$$x_q$$

$$y_h$$

$$K$$

$$=\frac{\left(1-p_A\right)\left(1-p_B\right)}{p_A}\sum _{h=1}^j{g}_h{B}_h.$$

$$\sum _{q=1}^i{f}_q{m}_q=\sum _{q=1}^i{m}_q+\sum _{h=1}^j{n}_h-1$$

$$\sum _{h=1}^j{g}_h{n}_h=\sum _{h=1}^j{n}_h+\sum _{q=1}^i{m}_q-1.$$

$${\left(1-p_A\right)}^{{\sum }_{q=1}^i{f}_q{m}_q-{\sum }_{q=1}^i{m}_q-{\sum }_{h=1}^j{n}_h+1}$$

$${\left(1-p_B\right)}^{{\sum }_{h=1}^j{g}_h{n}_h-{\sum }_{h=1}^j{n}_h-{\sum }_{q=1}^i{m}_q+1}.$$

$$N\{m,n\}$$

$$\times \left(\sum _{h=1}^j{g}_h{n}_h-\sum _{h=1}^j{n}_h\right)!\prod _{q=1}^i\frac{x_q^{m_q}}{m_q!}\prod _{h=1}^j\frac{y_h^{n_h}}{n_h!}$$

$$x_q$$

$$y_h$$

$$K$$

$$\log n!\approx \log \left(\sqrt{2\pi n}\right)+n\log \left(\frac{n}{e}\right)+\log \left(1+\frac{1}{12n}\right).$$

$$\log N\{m,n\}$$

$$+\log \left(\sum _{h=1}^j{g}_h{n}_h-\sum _{h=1}^j{n}_h\right)!$$

$$+\sum _{q=1}^i\left(m_q\log x_q-\log m_q!\right)$$

$$+\sum _{h=1}^j\left(n_h\log y_h-\log n_h!\right).$$

$$m_1+2m_2$$

$$m_2=n_1+n_2-1.$$

$$m_1=n_2+2.$$

$$P_l(L_l+R_l\to L_l{R}_l)=k_l{n}_{L_l}{n}_{R_l},$$

$$\begin{array}{rl}& n_{L_1}=n_{L_2}=2n_{\text{BPADA}},\\ & n_{L_3}=n_{L_4}=n_{\text{PA}},\\ & n_{R_1}=n_{R_3}=2n_{\text{MPD}},\\ & n_{R_2}=n_{R_4}=3n_{\text{TAPE}},\end{array}$$

$$\begin{array}{rl}& n_{L_1}=n_{L_2}=n_{\text{BPADA}},\\ & n_{L_3}=n_{L_4}=n_{\text{PA}},\\ & n_{R_1}=n_{R_3}=n_{\text{MPD}},\\ & n_{R_2}=n_{R_4}=n_{\text{TAPE}}.\end{array}$$

$$\alpha$$

$$=p_A{p}_B\frac{(1-{\rho }_1){\rho }_2}{1-p_A{p}_B(1-{\rho }_1)(1-{\rho }_2)},$$

$$\alpha =\frac{(1-{\rho }_1){\rho }_2}{{\rho }_1+{\rho }_2-{\rho }_1{\rho }_2}.$$

Full text

I INTRODUCTION

The molecular weight distribution and architecture are two important characteristics of a system of polymer chains Rubinstein and Colby (2003). They strongly affect material properties such as dynamic moduli, fracture toughness, glass transition temperature, and viscosity Rubinstein and Colby (2003); Nunes et al. (1982); Suneel et al. (2002). Experimental methods for an accurate determination of molecular weight distributions are thus of great interest Suneel et al. (2002); Mead (1994); Williamson et al. (2016). Theoretically, it is also highly desirable if the molecular weight distribution of a polymer can be predicted a priori based on the knowledge of the polymerization reaction without even synthesizing the polymer. Such a theoretical method will be a valuable tool not only useful for understanding experimental measurements but also beneficial for other theories and models aiming to predict polymer properties. For example, Nichetti and Manas-Zloczower proposed a theoretical model to predict viscosity of a polymer melt based on its molecular weight distribution that was determined by fitting the gel permeation chromatography data to statistical distribution functions Nichetti and Manas-Zloczower (1998). A method to quickly generate molecular weight distributions before polymers are synthesized thus might be able to advance the predictive capability of such theories.

A theory on the constitution and molecular size distribution of a step-growth polymer was proposed by Flory and Stockmayer many years ago Flory (1941a, b, c); Stockmayer (1943, 1944, 1952) and has been frequently used to determine the gel point. Flory studied the polymerization of bifunctional monomers mixed with trifunctional and tetra functional branching units and made two fundamental assumptions. First, the same functional group has the same probability to react with another group and this probability is not affected by the length of the polymer to which the functional group belongs as well as the position of the functional group on that polymer. Secondly, ring polymers are not formed. Stockmayer extended the theory to branching units with arbitrary functionalities and derived the Stockmayer formula for the number of a polymer with a given composition, though ring structures were still excluded. The predictions of the Flory-Stockmayer theory, including the gel point and the average molecular weight, have been tested experimentally Peebles (1971); Matsumoto (1995); Matsumoto et al. (1999); Bannister et al. (2006); Gao et al. (2007); Schultz et al. (2009); Rosselgong et al. (2009). However, the entire distribution is hard to probe experimentally and often only some average molecular weight is measured. Practically, it is also difficult to directly predict the molecular weight distribution using the Flory-Stockmayer theory because of mathematical complexity of computing the numbers of all possible molecules in a branched polymer system. Furthermore, the Flory-Stockmayer theory is only expected to be valid below the gel point. Beyond the gel point, the formation of cyclics and closed loops in network structures (i.e., branching) becomes significant and the theory may fail Lyu et al. (2018a, b).

Monte Carlo (MC) methods are a class of techniques based on random sampling to numerically solve problems that have a probabilistic interpretation Landau and Binder (2005). MC methods have broad applications in polymer science Hsu and Grassberger (2011); Brandão et al. (2015), especially in polymer reaction engineering Brandão et al. (2015). Johnson and O’Driscoll used MC simulation to study sequence distributions in step-growth copolymerization Johnson and O’Driscoll (1984). Tobita applied MC simulation to a wide range of polymerization problems, including free-radical cross-linking copolymerization Tobita (1993), emulsion polymerization Tobita (1995a), the modification of polymer via crosslinking and degradation Tobita (1995b), long-chain branching and random scission Tobita (2001), and living radical polymerization Tobita (2006); Tobita and Yanase (2007). Hadicke and Stutz used an amine-cured epoxy as an example to compare the structure of step-growth networks obtained by MC simulation to that predicted by a branching theory Hädicke and Stutz (2002). He et al. applied a MC method to simulate self-condensing vinyl polymerization in the presence of multifunctional initiators and probed the role of reactive rate constants He et al. (2001, 2003). Rouault and Milchev Rouault and Milchev (1997) and He et al. He et al. (1997) performed MC simulations to study the kinetics and chain length distributions in living polymerization. Prescott used a MC model to show that chain-length dependent termination plays a significant role in living/controlled free-radical polymerization systems containing reversible transfer agents Prescott (2003). In a series of papers, Al-Harthi et al. used dynamic MC simulations to study atom-transfer radical polymerizations Al-Harthi et al. (2006, 2007); Al-Harthi et al. (2009a, b). Polanowski et al. Polanowski et al. (2010, 2011) and Bannister et al. Bannister et al. (2009) used MC methods to study the branching and gelation in living copolymerization of vinyl and divinyl monomers. Recently, Lyu et al. used a similar model to study the atom transfer radical and the conventional free radical polymerization of divinyl monomers and checked the applicability of the Flory-Stockmayer theory in such systems Lyu et al. (2018a, b). Gao et al. used kinetic MC methods to simulate free radical copolymerization processes and discussed how to accelerate such simulations using scaling approaches Gao et al. (2015, 2017). Meimaroglou et al. proposed a MC algorithm to calculate the molecular weight distribution for linear polymers and the bivariate molecular weight-long chain branching distribution for highly branched polymers Meimaroglou et al. (2007). They also used MC simulation to investigate the molecular, topological, and solution properties of highly branched low-density polyethylene Meimaroglou et al. (2011) and the ring-opening homopolymerization of L,L-Lactide Meimaroglou et al. (2017). Iedema et al. developed a MC simulation model including both branching and random scission to calculate the molecular weight and branching distributions and compared their calculations to experimental measurements on high-density polyethylene Iedema et al. (2013). Yaghini and Iedema compared the results on low-density polyethylene from such MC simulations to the predictions of a multiradical model based on a Galerkin finite element approach Yaghini and Iedema (2015).

One important application of MC simulations is to quickly compute molecular weight distributions Tobita (1993, 1995a, 1995b, 2001, 2006); Tobita and Yanase (2007); Rouault and Milchev (1997); Milchev et al. (2000); He et al. (1997); Gao and He (2015); Maafa et al. (2007); Soares and Hamielec (2007); Hsu and Grassberger (2011); Hamzehlou et al. (2013); Iedema et al. (2013); Yaghini and Iedema (2015); Gao et al. (2015, 2017); Lyu et al. (2018a, b). In MC simulations, all structures including rings and networks allowed by a polymerization reaction can be produced Fawcett et al. (1995), no matter the system is below or beyond the gel point Tripathi et al. (2015). Various schemes can also be implemented for the kinetics of the polymerization, which thus allows us to test the specific assumptions made by a theory. The Gillespie algorithm can be used to speed up the kinetics of a reaction in silico and enable a reactive system to quickly reach a steady state Gillespie (1977, 2007). In this paper, we develop a MC simulation model based on the Gillespie algorithm to study the polymerization of polyetherimides (PEIs) in the presence of chain terminators and branching agents. The results from the MC simulations are used to test the Flory-Stockmayer theory including its assumption on the reaction rates.

This paper is organized as follows. In Sec. II the Flory-Stockmayer theory is introduced, the technical challenge of computing molecular weight distributions with this theory is discussed, and an approximation method is proposed. In Sec. III we describe the MC model of the polymerization process of branched PEIs in detail. In Sec. IV, MC results are compared to the predictions of the Flory-Stockmayer theory. Practically, computations of molecular weight distributions can be only be executed for a small system either with the Flory-Stockmayer theory or the MC model. We thus include a discussion on the effect of finite system size in this section. Although the emphasis is on stoichiometric, fully reacted systems, those that are only partially reacted and/or nonstoichiometric are also discussed in Sec. IV. Finally, a brief summary is provided in Sec. V.

II Flory-Stockmayer Theory of Step-Growth Polymers

Flory and Stockmayer considered a general step-growth polymer that consists of two types of monomers, $A$ and $B$. All reactions occur between $A$ and $B$. There are $i$ type-$A$ monomers denoted as $A_{1}$, $A_{2}$, …, $A_{i}$. To simplify notation, we also use $A_{q}$ with $q\in\{1,2,...,i\}$ to denote the number of $A_{q}$ monomers. Similarly, there are $j$ type-$B$ monomers and the corresponding numbers are $B_{1}$, $B_{2}$, …, $B_{j}$, respectively. The symbol $f_{q}$ denotes the functionality of an $A_{q}$ monomer, where $q\in\{1,2,...,i\}$, i.e., there are $f_{q}$ functional groups on an $A_{q}$ monomer that can form bonds with the corresponding functional groups on a $B_{h}$ monomer, where $h\in\{1,2,...,j\}$. The functionality of a $B_{h}$ monomer is denoted as $g_{h}$. The Flory-Stockmayer theory can be applied to a polymerized state where a fraction $p_{A}$ of all the functional groups on the type-$A$ monomers have reacted with a fraction $p_{B}$ of all the functional groups on the type-$B$ monomers. Therefore,

[TABLE]

In this paper, we call the systems with $\sum_{q=1}^{i}f_{q}A_{q}=\sum_{h=1}^{j}g_{h}B_{h}$ and thus $p_{A}=p_{B}$ as stoichiometric systems while those with $\sum_{q=1}^{i}f_{q}A_{q}\neq\sum_{h=1}^{j}g_{h}B_{h}$ and $p_{A}\neq p_{B}$ as nonstoichiometric. The systems with $p_{A}$ or $p_{B}$, or both, equal to 1 are fully reacted.

We use $N\{m,n\}$ to denote the number of molecules formed by $m_{q}$ monomers of sub-type $A_{q}$ and $n_{h}$ monomers of sub-type $B_{h}$, with $q$ running from 1 to $i$ and $h$ running from 1 to $j$. Here $\{m,n\}$ is a shorthand of $\{m_{1},m_{2}.,...,m_{i},n_{1},n_{2},...,n_{j}\}$, which denotes the monomer composition of a given molecule. The Flory-Stockmayer theory predicts that

[TABLE]

with

[TABLE]

Equation (2) is called the Stockmayer formula, which gives the number of molecules of any monomer compositions. However, practically it is difficult to compute the molecular weight distribution from the Stockmayer formula, as all the possible combinations for $\{m_{1},m_{2}.,...,m_{i},n_{1},n_{2},...,n_{j}\}$ have to be taken into account. Since $m_{q}$ runs from 1 to $A_{q}$ for $q\in\{1,2,...,i\}$ and $n_{h}$ runs from 1 to $B_{h}$ for $h\in\{1,2,...,j\}$, the total number of possible molecules is $\prod_{q=1}^{i}A_{q}!\times\prod_{h=1}^{j}B_{h}!$. This number is huge when there are many sub-types (i.e., large $i$ and $j$) and/or large numbers (i.e., large $A_{q}$ and $B_{h}$) of monomers involved in a polymerization.

For a molecule with composition $\{m,n\}$, the total number of monomers is $\sum_{q=1}^{i}m_{q}+\sum_{h=1}^{j}n_{h}$. Since the Flory-Stockmayer theory does not consider rings, the total number of bonds in this molecule must be $\sum_{q=1}^{i}m_{q}+\sum_{h=1}^{j}n_{h}-1$. When $p_{A}=p_{B}=1$, all the functional groups have reacted and in a given molecule the total number of the functional groups on all the type-$A$ monomers is equal to the total number of the functional groups on all the type-$B$ monomers. This number must also be equal to the total number of bonds in that molecule. Namely, for $p_{A}=p_{B}=1$ there are two identities,

[TABLE]

and

[TABLE]

These two identities can help us simplify the Stockmayer formula for stoichiometric, fully reacted systems. Note that the terms in Eq. (2) involving $\left(1-p_{A}\right)$ and $\left(1-p_{B}\right)$ appear as

[TABLE]

and

[TABLE]

These terms can be dropped out because of Eqs. (6) and (7). As a result, for fully reacted stoichiometric systems with $p_{A}=p_{B}=1$ the Stockmayer formula is simplified as

[TABLE]

with

[TABLE]

Computing $N\{m,n\}$ is not easy as it contains many factorials. The calculation can be expedited using the Stirling approximation,

[TABLE]

Then for fully reacted stoichiometric systems, the Stockmayer formula can be approximated logarithmically as

[TABLE]

The computation of the molecular weight distribution from $N\{m,n\}$ can be further accelerated by noting that not all the combinations $\{m,n\}$ will yield a molecule. For a fully reacted stoichiometric system where $p_{A}=p_{B}=1$, Eqs. (1), (6), and (7) can be used to reduce the total number of $\{m,n\}$. For the branched PEIs considered in this paper (see Sec. III), $f_{1}=1$, $f_{2}=2$, $g_{1}=2$ and $g_{2}=3$. The constraints become

[TABLE]

and

[TABLE]

Combining Eqs. (14) and (15), we get

[TABLE]

Equations (15) and (16) indicate that $m_{1}$ and $m_{2}$ are totally constrained by $n_{1}$ and $n_{2}$ in an allowed composition. Furthermore, since $m_{2}\geq 0$, $n_{1}$ and $n_{2}$ cannot be zero at the same time. The time complexity to enumerate all possible molecules is thus $\mathcal{O}(B_{1}B_{2})$, which is about $\mathcal{O}(Z^{2})$ with $Z$ being the system size (i.e., the total number of monomers prior to polymerization). This time complexity is acceptable for small systems. However, if there are more sub-types of monomers, then the time complexity will increase exponentially as $\mathcal{O}(Z^{w})$, where $w$ is the number of monomer sub-types. For not fully reacted or nonstoichiometric systems where $p_{A}$ or $p_{B}$ are less than 1, we lose the constrains that help reduce the number of possible $\{m,n\}$ and then computing the molecular weight distribution from $N\{m,n\}$ has to rely on Eq. (2) and will become more challenging, even though the Stirling approximation may still be used. In these situations, the MC model described below will serve as a solution as it does not suffer from such limitations and the time complexity of computing the molecular weight distribution with MC simulations is always $\mathcal{O}(Z)\times k$, where $k$ is the number of MC runs needed to obtain desired statistics. Usually, $k$ is about $10^{3}$ to $10^{4}$.

III Monte Carlo Model of Polymerization of Branched Polyetherimides

Four types of monomers are involved in the formation of branched PEIs, including 4,4’-bisphenol A dianhydride (BPADA), m-phenylenediamine (MPD), phthalic anhydride (PA), and tris[4-(4-aminophenoxy)phenyl] ethane (TAPE).Odle et al. The chemical structures of these monomers are shown in Fig. 1. The involved reaction is the condensation reaction between an amine group on MPD or TAPE and a carboxylic anhydride group on BPADA or PA. In the notation of the Flory-Stockmayer theory, PA is monomer $A_{1}$ with $f_{1}=1$, BPADA is monomer $A_{2}$ with $f_{2}=2$, MPD is monomer $B_{1}$ with $g_{1}=2$, and TAPE is monomer $B_{2}$ with $g_{2}=3$. Out of these monomers, PA is an end capper to terminate a chain and TAPE is a trifunctional branching agent. Fig. 1 shows the representation of these monomers in our MC model. Each functional group containing one carboxylic anhydride is mapped to an $A$ bead and that containing one amine is mapped to a $B$ bead. Each $A$ bead can react with a $B$ bead to form a bond (i.e., $A+B\rightarrow AB$), which describes the condensation reaction between an amine group and a carboxylic anhydride group.

In the formation of branched PEIs consisting of the above 4 types of monomers, there are 4 possible reactions, as sketched in Fig. 2. Reaction 1 is between BPADA and MPD, which leads to the formation of a PEI backbone. Reaction 2 is between BPADA and TAPE that results in branching. Reaction 3 is between PA and MPD, which terminates a PEI chain. Reaction 4 is between PA and TAPE, which consumes one amine group on TAPE and effectively reduces its functionality by 1.

With the mapping in Fig. 1 and the reaction scheme in Fig. 2, we perform MC simulations to study the polymerization of branched PEIs. We adopt the Gillespie algorithm to speed up MC simulations. Since only the final chain structures are concerned, we neglect the random process in the typical Gillespie algorithm that determines the time interval after which the next reaction occurs. We only keep the random process of picking a reaction at a time. At each MC step, all 4 reactions will have a probability to occur and the reaction rate of a particular reaction is determined by a rate constant and the concentration of the two types of monomers involved in that reaction. Mathematically, the probability of reaction $l$ is proportional to

[TABLE]

where $L_{l}$ ($R_{l}$) represents the reactant consisting of $A$ ($B$) beads, $k_{l}$ is a rate constant, $n_{L_{l}}$ ($n_{R_{l}}$) is a quantity that depends on the concentration of the reactant $L_{l}$ ($R_{l}$), and $l\in\{1,2,3,4\}$ indexes the possible reactions sketched in Fig. 2. Specifically, $L_{1}$ and $L_{2}$ are BPADA, $L_{3}$ and $L_{4}$ are PA, $R_{1}$ and $R_{3}$ are MPD, and $R_{2}$ and $R_{4}$ are TAPE. Since all the 4 reactions can be reduced to the reaction between an $A$ bead and a $B$ bead (i.e., the reaction between a functional group containing one amine and another functional group containing one carboxylic anhydride) as shown in Fig. 1(e), $k_{l}$ will be set as a constant $k$ for all the 4 reactions and $n_{L_{l}}$ and $n_{R_{l}}$ can be expressed as

[TABLE]

where $n_{\text{BPADA}}$, $n_{\text{PA}}$, $n_{\text{MPD}}$, and $n_{\text{TAPE}}$ are the concentrations of monomers available for reactions (i.e., monomers with at least one unreacted functional group). In other words, $n_{L_{l}}$ ($n_{R_{l}}$) is the concentration in terms of the number of $A$ ($B$) beads on the reactant $L_{l}$ ($R_{l}$). The particular reason of writing $n_{L_{l}}$ and $n_{R_{l}}$ in this way will be discussed in Sec. IV.1.

At each MC step, the probability of Reaction $l$ to be chosen is equal to $P_{l}/\sum_{q=1}^{4}P_{q}$. After a reaction is selected, a pair of $L_{l}$ and $R_{l}$ that have unreacted functional groups (i.e., with unreacted $A$ and $B$ beads, respectively) is randomly chosen to react. Then the system status is updated, including the bond information between the monomers and the identity of monomers with unreacted functional groups. The MC process is repeated for the updated system until no more reactions can occur or the system has reached a desired extent of reaction. The flow chart of the MC simulation model is shown in Fig. 3. Note that in this model, we made a simplification by not allowing backward reactions, which means that once formed, the bond between an $A$ bead and a $B$ bead cannot break. However, the model permits the formation of both rings and networks.

IV Results and Discussion

IV.1 Rate Constant $k$

Equation (18) indicates that the reaction rate $P_{l}$ is based on the concentrations of the functional groups (i.e., $A$ beads or $B$ beads) on the reactants involved in that reaction. However, $P_{l}$ can also be computed from the concentrations of the available reactants themselves, i.e., the monomer concentrations. In this case, the reaction rate $P_{l}$ can be written in the same way as in Eq. (17) but with Eq. (18) replaced by

[TABLE]

To check which way of computing the reaction rates yields results that are more applicable to realistic systems, we performed a test with a simple system consisting of only PA and TAPE monomers, as shown in Table 1. For this system, there are only 4 possible final products, including single TAPEs and TAPEs connected with one, two, or three PAs, respectively.

Figure 4 shows the results on the probability distribution of the 4 final products for the system in Table 1, for which gelation is not a concern. The comparison shows that the results from the MC simulations based on Eq. (18) agree with the Flory-Stockmayer theory while those based on Eq. (19) do not. Furthermore, for the molar ratio in Table 1, all the anhydride groups on the PA monomers are reacted and each amine group on a TAPE monomer has a 2/3 chance to be reacted in a fully reacted system. A simple statistical analysis shows that the probabilities for a TAPE monomer to react with 0, 1, 2, and 3 PAs are $\left(\frac{1}{3}\right)^{3}$, $3\times\frac{2}{3}\times\left(\frac{1}{3}\right)^{2}$, $3\times\left(\frac{2}{3}\right)^{2}\times\frac{1}{3}$, and $\left(\frac{2}{3}\right)^{3}$, respectively. These results are also plotted in Fig. 4 and are close to those from the Flory-Stockmayer theory and the MC simulations based on Eq. (18). The small differences are due to the fact that the theory and simulations consider a finite system while the statistical model assumes an infinite system. We conclude that the reactions rates based on Eq. (18) should be used in the MC simulations. From now on all the data presented in this paper are computed with Eq. (18) for the reaction rates. In the next two sections, Sec. IV.2 and Sec. IV.3, we focus on fully reacted stoichiometric systems where $p_{A}=p_{B}=1$. We discuss partially reacted stoichiometric systems where $p_{A}=p_{B}<1$ in Sec. IV.4 and nonstoichiometric systems where $p_{A}\neq p_{B}$ in Sec. IV.5.

IV.2 Fully Reacted Stoichiometric Systems

For the branched PEIs considered in this paper, type $A$ monomers are BPADA and PA and type $B$ monomers are MPD and TAPE, with $f_{1}=1$, $f_{2}=2$, $g_{1}=2$ and $g_{2}=3$. From the Flory theory Flory (1941a), the gel point is $\alpha_{c}=1/(g_{2}-1)=1/2$. However, for the systems at hand the expression of $\alpha$, which characterizes the level of cross-linking, has to be modified because each PA monomer as a chain terminator has only one functional group. The modified expression is

[TABLE]

where $\rho_{1}$ is the fraction of functional groups on the terminators (i.e., PA monomers) with respect to all the functional groups on type $A$ monomers and $\rho_{2}$ is the fraction of functional groups on the branching agents (i.e., TAPE monomers) with respect to all the functional groups on type $B$ monomers. For a fully reacted stoichiometric system, $p_{A}=p_{B}=1$ and this expression can be simplified as

[TABLE]

We can vary the numbers of monomers to tune $\rho_{1}$ and $\rho_{2}$, thus putting the fully reacted system below, around, or beyond the gel point. Three such systems are listed in Table 2, where $\rho_{2}$ is changed by varying the numbers of MPD and TAPE monomers. In the MC simulations of these stoichiometric systems, we set $p_{A}=p_{B}=1$, thus allowing the systems to be fully reacted.

The results on the molecular weight distribution from the Flory-Stockmayer theory and the MC simulations are shown in Fig. 5. The comparison shows that for a system below the gel point such as $S_{<}$ (Fig. 5(a)), the MC results agree well with the Flory-Stockmayer theory. For $S_{\simeq}$ which is close to the gel point, the Flory-Stockmayer theory overestimates the fraction of low molecular weight polymers and underestimates the fraction of high molecular weight species when compared to the results from the MC simulations, as shown in Fig. 5(b). The discrepancy between the Flory-Stockmayer theory and the MC results becomes more dramatic for systems above the gel point. For $S_{>}$, $\alpha=0.74308$, way above the critical gel point $\alpha_{c}=0.5$. The Flory-Stockmayer theory predicts a probability density that is about 8 times of the MC result in the region of low molecular weight from 0 to about $0.5\times 10^{5}$ Da, as shown in Fig. 5(c). However, the MC simulations show a significant fraction of polymers in the region of molecular weight higher than about $1.5\times 10^{5}$ Da and these high molecular weight polymers are completely overlooked by the Flory-Stockmayer theory, as shown in the inset of Fig. 5(c). This discrepancy is not surprising as beyond the gel point, polymers with a large network structure are expected and closed loops can frequently emerge in such polymers. The Flory-Stockmayer theory does not consider the formation of rings and thus cannot accurately predict the molecular weight distribution for systems above the gel point.

IV.3 Effect of System Size

In experiments, the amount of monomers involved is at the order of moles, i.e., at the order of $10^{23}$. It is thus practically impossible to directly compute the molecular weight distribution from the Stockmayer formula (Eq. (2)) for such macroscopic systems. These systems are also out of the reach of MC simulations that typically deals with systems of fewer than $10^{6}$ monomers. A natural question we can ask is: if we keep the molar ratios unchanged but reduce the numbers of participating monomers in proportion, can we use either the Flory-Stockmayer theory or the MC simulations to generate a molecular weight distribution that is applicable to a macroscopic system? To answer this question, we test 4 additional systems listed in Table 3. The smallest system has 10 PA, 134 BPADA, 146 MPD, and 2 TAPE and is denoted as $S_{1}$. Then the numbers of monomers are increased 10, 50, and 80 fold by keeping the ratios to generate systems $S_{10}$, $S_{50}$, and $S_{80}$. The subscript of the system label reflects the size ratio with respect to the smallest system, $S_{1}$. In this notation, the system $S_{<}$ in Table 2 is equivalent to $S_{5}$. All these systems are still below the gel point when fully reacted.

The molecular weight distributions predicted by the Flory-Stockmayer theory for $S_{1}$, $S_{<}$, $S_{10}$, $S_{50}$, and $S_{80}$, including the probability density and the cumulative probability, are shown in Fig. 6. The main panels are for the region of low molecular weight and the insets show the data in the high molecular weight region. The data show that when the system size is increased, the curves of the molecular weight distribution converge quickly. There is a clear difference between the data for $S_{1}$ and those for $S_{<}$ (i.e., $S_{5}$). However, the difference between $S_{<}$ and $S_{80}$ is very small in the low molecular weight region and only discernible in the tail of the distribution in the region of high molecular weight (see the insets of Fig. 6). Furthermore, the results for $S_{50}$ and $S_{80}$ are almost indistinguishable in the entire region of molecular weight relevant to experiments, indicating that these systems are already large enough such that the molecular weight distribution is not affected by the finite system size any more.

Since $S_{1}$, $S_{<}$, $S_{10}$, $S_{50}$, and $S_{80}$ are all below the gel point, we expect the results from the Flory-Stockmayer theory and the MC simulations on the molecular weight distribution to agree. The comparison between the two is shown in Fig. 7 for $S_{1}$ and $S_{80}$. For $S_{1}$, some difference is observed between the prediction of the Flory-Stockmayer theory and the MC result because of the small size of this system. An excellent agreement is found between the theory and simulations for $S_{80}$. Similar agreements are also found for $S_{10}$ and $S_{50}$. A good agreement is already discussed earlier for $S_{<}$ as shown in Fig. 5(a). These comparisons once again confirm that the Flory-Stockmayer theory provides a good description of the molecular weight distribution for systems well below the gel point, where ring formation is not a big concern. Below the gel point, both the Flory-Stockmayer theory and the MC simulations can be applied to a system containing only a few hundred to a few thousand monomers but having the same molar ratios of monomers as a macroscopic system to accurately predict the molecular weight distribution. As discussed earlier, the Flory-Stockmayer theory starts to fail when a system approaches or goes above the gel point. However, in these situations the MC simulations can still be used to quickly generate a molecular weight distribution that is applicable to an experimental system.

IV.4 Partially Reacted Stoichiometric Systems

Up to this point, we mainly focus on fully reacted stoichiometric systems as it is possible to compute the molecular weight distribution using the Stockmayer formula even for a system with a relatively large size such as $S_{80}$. In this and next section we show that the conclusions reached so far also apply to partially reacted and/or nonstoichiometric systems. However, because of the practical difficulty of using the Stockmayer formula to compute the molecular weight distribution when either $p_{A}$ or $p_{B}$, or both, are less than 1, we use small systems with sizes similar to $S_{5}$ to illustrate the main point.

In this section we discuss partially reacted stoichiometric systems where $\sum_{q=1}^{i}f_{q}A_{q}=\sum_{h=1}^{j}g_{h}B_{h}$ but $p_{A}=p_{B}<1$. Five such systems with the same size as $S_{>}$ are listed in Table 4 where the values of $p_{A}$ and $p_{B}$ are increased from 0.95 to 0.99. The corresponding values of $\alpha$ changes from about 0.42 to about 0.65, thus enclosing the gelation transition at $\alpha_{c}=0.5$.

The results on the molecular weight distribution from the Flory-Stockmayer theory and MC simulations at various values of $p_{A}$ and $p_{B}$ are shown in Figs. 8(a) and (b), respectively. The molecular weight distribution predicted by the Flory-Stockmayer theory seems to be relatively insensitive to the values of $p_{A}$ and $p_{B}$. However, the MC results show that when the value of $p_{A}$ and $p_{B}$ is increased, the probability density in the low molecular weight region is reduced (see Fig. 8(b)) while that in the high molecular weight region is enhanced (see the inset of Fig. 8(b)). This systematic trend is expected as when the extent of reaction is larger, more polymers with higher molecular weights are anticipated to form.

To compare the predictions of the Flory-Stockmayer theory to the MC results on the molecular weight distribution, in Fig. 8(c) their differences are shown for various $p_{A}$ and $p_{B}$. It is clear that when $p_{A}$ and $p_{B}$ are small, the systems are below the gel point and the results from the theory and simulations agree, as for $S^{0.95}$ and $S^{0.96}$. The difference becomes noticeable when the system approaches the gel point, such as $S^{0.97}$. For $S^{0.98}$ and $S^{0.99}$, they are above the gel point and clear differences in the probability density from the theory and simulations can be noted in both low (see Fig. 8(c)) and high (see the inset of Fig. 8(c)) molecular weight regions. The results for the partially reacted stoichiometric systems thus corroborate the conclusion that the Flory-Stockmayer theory only applies to systems well below the gel point. However, the MC simulations can be used to compute the molecular weight distribution for any systems no matter they are below, around, or above the gel point.

IV.5 Nonstoichiometric Systems

We finally discuss nonstoichiometric systems where $\sum_{q=1}^{i}f_{q}A_{q}\neq\sum_{h=1}^{j}g_{h}B_{h}$ and $p_{A}\neq p_{B}$. Three systems with sizes similar to $S_{>}$ are shown in Table 5. We fix the value of $p_{A}$ at 0.99 but vary $p_{B}$ from 0.93 to 0.97. The numbers of monomers of PA, BPADA, and TAPE are all fixed. The number of MPD is varied according to Eq. (1). Specifically, when the number of MPD is reduced, the values of $p_{B}$, $\rho_{2}$, and $\alpha$ are all increased. For the three systems in Table 5, $S_{<}^{n}$ is below, $S_{\simeq}^{n}$ is around, and $S_{>}^{n}$ is above the gel point. Here the superscript $n$ in the system labels indicates that these systems are nonstoichiometric.

The results on the molecular weight distribution for the three nonstoichiometric systems are plotted in Fig. 9. For $S_{<}^{n}$ which is below the gel point, the MC results agree with the prediction of the Flory-Stockmayer theory, as shown in Fig. 9(a). The two start to differ when a system approaches the gel point. An example is shown in Fig. 9(b) for $S_{\simeq}^{n}$ with $\alpha=0.502$. For this system, the Flory-Stockmayer theory overestimates the probability of low molecular weight polymers while underestimates the probability in the region of molecular weight higher than about $0.5\times 10^{5}$ Da (see the inset of Fig. 9(b)). For $S_{>}^{n}$ which is above the gel point, the MC results on the probability density are smaller than those calculated with the Flory-Stockmayer theory when the molecular weight is lower than about $1\times 10^{5}$ Da (Fig. 9(c)) but higher at higher molecular weights (see the inset of Fig. 9(c)). For $S_{>}^{n}$ the molecular weight distribution has a second peak around $2\times 10^{5}$ Da, while the Flory-Stockmayer theory predicts a monotonically decaying distribution in this region. The results on nonstoichiometric systems thus one more time indicate that the Flory-Stockmayer theory only applies to systems well below the gel point, for which the formation of cyclic polymers or closed loops is negligible.

V Conclusions

We have used MC simulations to study the polymerization of branched PEIs from BPADA (backbone monomer), MPD (backbone monomer), PA (chain terminator), and TAPE (branching agent). All the reactions for this system can be reduced to a condensation reaction between an amine group and a carboxylic anhydride group and thus can be characterized by one reaction rate. Our work show that in the MC model, the reaction rate should be computed using the concentrations of the functional groups on the monomers involved in a specific reaction, not the concentrations of the monomers themselves. The MC results are compared to the predictions of the Flory-Stockmayer theory. A practical approach of using the Flory-Stockmayer theory to compute molecular weight distributions has been suggested. We find that both the Flory-Stockmayer theory and the MC simulations accurately predict the molecular weight distribution for systems well below the gel point that is set by the functionality of the branching agent, though ring formation is not considered by the Flory-Stockmayer theory but allowed in MC simulations. The agreement between the theory and simulations thus indicates that ring formation is negligible for systems well below the gel point. However, for systems close to or above the gel point, the Flory-Stockmayer theory is not applicable as many cyclic polymers can be produced and ring structures can form in highly branched networks. For these systems, the MC simulations can still be used to quickly compute the molecular weight distribution that can be used to describe experimental measurements including average molecular weights.

Our tests indicate that in the MC simulations, a system with only a few hundred to a few thousand monomers but the same molar ratios of participating monomers is large enough to yield converging results on the molecular weight distribution for the region of molecular weight relevant to typical experiments (from 0 to about $3\times 10^{5}$ Da in the case of PEIs). These conclusions have been thoroughly confirmed with simulations for fully reacted, partially reacted, stoichiometric, and nonstoichiometric systems. The MC model presented here is expected to be applicable to a wide range of step-growth polymers.

Acknowledgement

This paper is based on the results from work supported by SABIC Innovative Plastics US LLC.