Deconvolution of GPC data

In my previous post, I developed a GPC post-processing program that uses the continuous wavelet transform (CWT) to detect and correct baselines, and the MH equation to convert GPC data into a molecular weight distribution. In this follow-up work, I determine the molecular weight distributions in a two-step synthesis process. In the first step, polymers are pre-polymerized with a distribution denoted as $P_1(M)$. In the second step, these polymers continue to grow, ultimately yielding a final chain-length distribution $Q(M)$. The goal here is to extract the molecular weight distribution for the second step, denoted as $P_2(M)$.

I present a straightforward method to calculate the molecular weight distribution of the second block using the GPC data from the preceding block and the final diblock copolymer. Assuming that we already know the distribution of the preceding block, $P_1(n)$, it is generally expected that the growth of the second block depends on the preceding chain length $n$. Consequently, the joint probability density function (pdf) for the two blocks is given by
$$ Q(n, m) = P_1(n)P_2(m, n)$$
and the pdf of the final diblock copolymer is obtained as
$$ Q(x) = \int Q(n, x-n) \mathrm{d}n $$

It is reasonable to assume that $P_2(m, n)$ can be factorized into the product of two functions, say $f(m)g(n)$, by neglecting higher-order correlations between $m$ and $n$. In this factorization, $g(n)$ acts as a scaling factor on the distribution $P_1(n)$, and the calculation of $P(x)$ becomes a convolution. The simplest case is when $g(n)=1$, meaning that the growth of the second block is independent of the length of the preceding block. Alternatively, one might assume $g(n) \sim n^{-1}$ to account for diffusion effects, implying that shorter preceding chains tend to grow a longer second block.

The evaluation of the pdf for the second block proceeds in three steps:

  1. Determine the range of molecular weights as $(x_{min} – n_{max},\, x_{max} – n_{min})$, where $x$ is the chain length of the final diblock copolymer and $n$ is the length of the preceding block.
  2. Interpolate the GPC-derived distributions onto an evenly spaced molecular weight grid, setting any negative values to zero.
  3. Deconvolute $Q(x)$ with $P_1(n)$.

This method provides a straightforward pathway to isolate the molecular weight distribution of the second synthetic block from the overall diblock copolymer distribution.

Phenomenological relation of coil to globule transition of polymer chain

Consider a polymer chain dissolved in a solvent, where its behavior is influenced by the Flory parameter, denoted by $\chi$. As the interaction parameter $\chi$ varies, the chain undergoes several conformational transitions:

• For $\chi > \chi_c$, the polymer collapses into a dense, globular structure with a radius of gyration scaling as $R_g \sim N^{1/3}$.

• At the critical point—commonly known as the $\theta$ point—where $\chi \sim \chi_c$, the chain behaves like an ideal random coil, characterized by $R_g \sim N^{1/2}$.

• In a good solvent, when $\chi < \chi_c$, the polymer chain is swollen and stretched, with $R_g \sim N^{3/5}$.

A “universal” expression that captures how the radius of gyration $R_g$ changes with $\chi$ is given by

  $R_g(\chi) = R_g^\mathrm{glob} + \frac{R_g^\mathrm{coil}-R_g^\mathrm{glob}}{1 + \exp[(\chi-\chi_c)/\Delta \chi]}$,

where $\Delta \chi$ quantifies the width of the conformational transition region and $\chi_c$ locates the transition.

To estimate the window width $\Delta \chi$, we first define an order parameter that measures deviations of the polymer’s size from its value at the $\theta$ point. We write

  $R = R_0 (1 + m)$,

with $R_0 \sim N^{1/2}b$, where $b$ is the monomer size and $m$ represents the fractional deviation from the ideal size.

Next, we expand the free energy $F(m)$ in powers of $m$. Close to the $\theta$ point, the expansion takes the form

  $\frac{F(m)}{k_BT} \simeq C_1 N (\chi-\chi_c) m^2 + C_2 N m^4 + \cdots$,

where $C_1$ and $C_2$ are constants, and $k_BT$ is the thermal energy.

In the mean-field (infinite-system) limit, the phase transition is sharp. However, for finite systems, fluctuations smear the transition, leading to a rounded behavior characterized by a finite width $\Delta \chi$ defined by $|\chi-\chi_c|$.

At the edge of the transition, the quadratic and quartic terms in the free energy become comparable. Equating these contributions gives

  $C_1 N\,\Delta\chi\,m^2 \sim C_2 N\,m^4,$

which implies

  $m^2 \sim \frac{C_1}{C_2}\,\Delta\chi$.

Furthermore, the transition is significantly broadened when the free energy barrier, estimated by the quartic term, is of the order of unity (in units of $k_BT$). That is, when

  $N C_2\, m^4 \sim O(1).$

Substituting our earlier estimate for $m^2$ into this condition, we obtain

  $N C_2 \left(\frac{C_1}{C_2}\Delta\chi\right)^2 \sim O(1),$

which simplifies to

  $\Delta\chi^2 \sim \frac{1}{N}.$

Thus, the rounding of the transition in terms of $\chi$ scales as

  $\Delta\chi \sim \frac{1}{\sqrt{N}}.$

Distribution of segments on Gaussian chain

This analysis focuses on the probability distribution function, $P_i(\mathbf{r}_i-\mathbf{r}_\mathrm{cm})$, of the $i$th segment relative to the center of mass in an ideal chain. An ideal chain is modeled as a multidimensional random walk with independent steps. The distribution of each step, with a mean length of $b$, is assumed to be Gaussian: $P(\mathbf{r})\sim\mathcal{N}(0,b^2)$.

Let $\mathbf{b}_i=\mathbf{r}_{i+1}-\mathbf{r}_{i}$ represent the $i$th bond vector. Then, the position of the $i$th segment is given by:

$\mathbf{r}_i=\sum_{j=1}^{i-1} \mathbf{b}_j$

The center of mass, $\mathbf{r}_\mathrm{cm}$, is calculated as:

$\mathbf{r}_\mathrm{cm}=\frac{1}{N}\sum_{i=1}^{N} \mathbf{r}_i = \frac{1}{N}\sum_{j=1}^{N-1}(N-j)\mathbf{b}_j$

Therefore, the displacement of the $i$th segment relative to the center of mass is:

$\mathbf{r}_i-\mathbf{r}_\mathrm{cm}=\sum_{j=1}^{i-1}\frac{j}{N}\mathbf{b}_j+\sum_{j=i}^{N-1}\frac{N-j}{N}\mathbf{b}_j$

If $X$ is a Gaussian random variable with variance $\sigma^2$, then $aX$ is also a Gaussian random variable with variance $a^2\sigma^2$. Using this property, we can write the characteristic function for $P_i(\mathbf{r}_i-\mathbf{r}_\mathrm{cm})$ of a $d$-dimensional ideal chain:

$\phi_i(\mathbf{q})=\Pi_{j} \phi_{\mathbf{b}_j’}(\mathbf{q})=\exp\left(-\frac{1}{2}\mathbf{q}^T\left(\sum_{j=1}^{N-1}\Sigma_j\right)\mathbf{q}\right)$

where $\mathbf{b}’_j=\frac{j}{N}\mathbf{b}_j$ for $j\le i-1$ and $\frac{N-j}{N}\mathbf{b}_j$ for $i\le j \le N-1$. $\phi_{\mathbf{b}_j’}=-\exp(-0.5\mathbf{q}^T\Sigma\mathbf{q})$ is the characteristic function of the probability distribution of bond $j$. $\Sigma_j=\frac{j^2 b^2}{d N^2}\mathbb{I}_d$ for $j\le i-1$ and $\Sigma_j=\frac{(N-j)^2 b^2}{d N^2}\mathbb{I}_d$ for $i\le j \le N-1$, where $\mathbb{I}_d$ is the $d$-dimensional identity matrix. Setting $b=1$ for convenience, we obtain:

$\phi_i(\mathbf{q})=$ $\exp \left(-\frac{\left(6 i^2-6 i (N+1)+2 N^2+3 N+1\right) \left(q_x^2+q_y^2+q_z^2\right)}{36 N}\right)$

The corresponding distribution of this characteristic function is still a Gaussian distribution with $\Sigma=\frac{b^2}{3} \mathbb{I}_3$, where the equivalent bond length $b^2=\frac{\left(6 i^2-6 i (N+1)+2 N^2+3 N+1\right)}{6 N}$. The 6th moment is calculated as $\frac{1}{N}\sum_{i=1}^N \langle(\mathbf{r}_i-\mathbf{r}_\mathrm{cm})^6\rangle=\frac{58 N^6-273 N^4+462 N^2-247}{1944 N^3}$. For large $N$, only the leading term, $\frac{29}{972} N^3$, is significant. For $N=20$, the result is $235.886$, which agrees with simulations. Another example is for $N=5$, where the $R_g^2$ is $0.8$ compared to the expected value of $5/6=0.8\dot{3}$.

Here is the simulation code:

ch = np.random.normal(size=(100000,20,3),scale=1/3**0.5)
ch[:,0,:]=0
ch = ch.cumsum(axis=1)
ch -= ch.mean(axis=1,keepdims=1)
m6 = np.mean(np.linalg.norm(ch,axis=-1)**6)

Saddle point method and end-to-end distribution of polymer models

Introduction

Saddle point methods are widely used to estimate integrals of the form:

$I = \int \exp\left(-f(x)\right)\mathrm{d}x$

where the function $f(x)$ can be approximated by the first two terms of its Taylor series expansion around a point $x_0$:

$f(x)\approx f(x_0) + f^\prime (x_0)(x-x_0) + \frac{1}{2}f^{\prime\prime}(x_0)(x-x_0)^2$

The integral is then approximated by its value at the saddle point, where $f^\prime (x_0)=0$ and $f^{\prime\prime}(x_0)>0$:

$\begin{align} I&\approx \int \exp\left(-f(x_0) – \frac{1}{2}f^{\prime\prime}(x_0)(x-x_0)^2\right) \mathrm{d}x\\ &=\exp(-f(x_0))\sqrt{\frac{2\pi}{f^{\prime\prime}(x_0)}} \end{align}$

Examples

  • Stirling’s Formula:

Using the definition of the Gamma function, we have:

$N!=\int_0^\infty \exp(-x)x^N\mathrm{d}x$

Let $f(x):=x-N\ln(x)$. For large $N$, the negative term is negligible. Solving $f^\prime (x) = 0$, we obtain:

$N!\approx\exp(-N+N\ln(N))\sqrt{2\pi{}N}=\sqrt{2\pi{}N}\left(\frac{N}{e}\right)^N$

  • Partition Function:

The partition function is given by:

$Z = \int \exp(-\beta U(\mathbf{x})) \mathrm{d}\mathbf{x}$

where $U(\mathbf{x})$ can be approximated by:

$U(\mathbf{x})\approx U(\mathbf{x}_0) + \frac{1}{2}(\mathbf{x}-\mathbf{x}_0)^T H[U](\mathbf{x}-\mathbf{x}_0)$

where $H$ represents the Hessian matrix.

End-to-End Distribution Function of Random Walk Model of Polymer Chains

For an $N$-step random walk model, the exact end-to-end vector distribution is:

$\begin{align} P(\mathbf{Y})&=\frac{1}{(2\pi)^3}\int \mathrm{d}\mathbf{k} \exp(-i\mathbf{k}\cdot\mathbf{Y})\tilde{\phi}^N\\ &=\int_{0}^{\infty} k\sin(kY) \left(\frac{\sin(kb)}{kb}\right)^N \mathrm{d}k \end{align}$

where $\phi(\mathbf{x})=\frac{1}{4\pi b^2}\delta(|\mathbf{x}|-b)$ is the distribution of a single step vector (length=$b$) and $\tilde{\phi}$ is its characteristic function; $\mathbf{Y}:=\sum_{i=1}^N \mathbf{x}_i$ is the end-to-end vector. Let $s=kb$ and $f(s):=i\frac{Y}{Nb}s -\ln\frac{\sin(s)}{s}$. Then we have:

$\begin{align} P&=\frac{i}{4\pi^2 b^2 Y}\int_{-\infty}^{+\infty} s \exp\left(-is\frac{Y}{b}\right)\left(\frac{\sin(s)}{s}\right)^N\mathrm{d}s\\ &=\frac{i}{4\pi^2 b^2 Y} \int s \exp(-Nf(s)) \mathrm{d} s\end{align}$

The integral is extended to $(-\infty, +\infty)$ due to the symmetry of the sine and cosine functions. The first sine function is replaced with an exponential form using Euler’s formula: $\exp(ix)=\cos(x)+i\sin(x)$.

Solving $f^\prime(s)=0$, we find that $is$ satisfies:

$\coth(is)-\frac{1}{is}=\frac{Y}{Nb}$

This is the Langevin function, $is_0=L^{-1}(\frac{Y}{Nb})$. Therefore, we have:

$\begin{align} P &\approx\frac{s_0}{4\pi^2b^2Y}\sqrt{\frac{2\pi}{Nf^{\prime\prime}(s_0)}}e^{-Nf(s_0)}\\ &=\frac{1}{(2\pi{}Nb^2)^{3/2}}\frac{L^{-1}(x)^2}{x(1-(L^{-1}(x)\csc\left(L^{-1}(x)\right))^2)^{1/2}}\\&\times\left(\frac{\sinh\left(L^{-1}(x)\right)}{L^{-1}(x)\exp\left(xL^{-1}(x)\right)}\right)^N\end{align}$

where $x:=\frac{Y}{Nb}$.

Steady compliance (Linear viscoelasty)

The Boltzmann superposition principle states that the strain response of a viscoelastic material is a superposition of the responses to all previous stress histories. Assuming no stress before time $t=0$, the constitutive equation can be written as:

$\begin{equation} \gamma(t)=\int_0^t J(t-t^\prime)\dot{\sigma}(t^\prime)\mathrm{d}t^\prime \end{equation}$

Applying the Laplace transform to this equation yields:

$\begin{align} \hat{\gamma}(s)&=\hat{J}(s)\hat{\dot{\sigma}}(s)\\ &=\hat{J}(s)\left(s\hat{\sigma}(s)-\sigma(0^{-})\right) \\ &=\hat{J}(s)\left(s\hat{G}(s)\hat{\dot{\gamma}}(s)-\sigma(0^{-})\right)\\ &=\hat{J}(s)\left(s\hat{G}(s)(s\hat{\gamma}(s)-\gamma(0^-))-\sigma(0^{-})\right)\\ &=s^2\hat{J}(s)\hat{\gamma}(s)\hat{G}(s)  \end{align}$

Here, $\hat{f}(s):=\mathcal{L}\lbrace f(t)\rbrace(s)$ denotes the Laplace transform of $f(t)$. Since the stress starts at time $t=0$, both $\gamma(0^-)$ and $\sigma(0^-)$ are zero. The second to third step is derived from the convolution relation $\sigma(t)=\int_0^t G(t-t^\prime)\dot{\gamma}(t^\prime)\mathrm{d}t^\prime$. Cancelling out $\hat{\gamma}(s)$ and rearranging the last equation gives:

$\begin{equation} \frac{1}{s^2}=\hat{J}(s)\hat{G}(s) \end{equation}$

Taking the inverse Laplace transform of this equation leads to the convolution relation:

$\begin{equation} t=\int_0^{t}J(t-t^\prime)G(t^\prime)\mathrm{d}t^\prime \end{equation}$

Substituting $J(t)=J_e + t/\eta$ into the convolution relation and taking the limit as $t\to\infty$ gives:

$\begin{equation} \begin{aligned} \color{blue}{t} &=\int_0^{t}\left(J_e+\frac{t-t^\prime}{\eta}\right)G(t^\prime)\mathrm{d}t^\prime\\ &= J_e\int_0^\infty G(t^\prime)\mathrm{d}t^\prime + {\color{blue}{\frac{t}{\eta}\int_0^\infty G(t^\prime)\mathrm{d}t^\prime}} – \int_0^\infty G(t^\prime)\frac{t^\prime}{\eta}\mathrm{d}t^\prime \end{aligned} \end{equation}$

Since $\eta:=\int_0^\infty G(t^\prime)\mathrm{d}t^\prime$, the blue terms cancel out, resulting in:

\begin{equation} J_e\int_0^\infty G(t^\prime)\mathrm{d}t^\prime = J_e\eta= \int_0^\infty G(t^\prime)\frac{t^\prime}{\eta}\mathrm{d}t^\prime = \frac{1}{\eta} \int_0^\infty G(t^\prime)t^\prime\mathrm{d}t^\prime \end{equation}

Therefore, we obtain:

$J_e = \frac{1}{\eta^2}\int_0^\infty tG(t)\mathrm{d}t$

$Q.E.D.$

Back2Top ^