This Blog is Untitled

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}$.

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.$

Persistence length, $L_p$, is a fundamental mechanical property that quantifies the bending stiffness of a polymer. It’s defined as the characteristic length scale over which the correlation of bond angles decays. This correlation is expressed as the average cosine of the angle, $\theta$, between bonds separated by a distance $s$ along the chain:

$\langle\cos(\theta(s))\rangle = \exp(-sl/L_p)$

To calculate $L_p$, we need to determine $\langle\cos(\theta(s))\rangle$. For adjacent bonds, where $s=1$, this simplifies to:

$\langle\cos(\theta)\rangle = \exp(-l/L_p)$

In simulations, the bending stiffness is often modeled using a bending energy, $U_b$. When $U_b$ is large, and excluded volume effects are negligible, we can calculate $\langle\cos(\theta)\rangle$ as:

$\langle\cos(\theta)\rangle=\frac{\int_0^\pi \cos(\theta)\sin(\theta)\exp(-\beta U_b)\mathrm{d}\theta}{\int_0^\pi\sin(\theta)\exp(-\beta U_b)\mathrm{d}\theta}$

The factor $\sin(\theta)$ represents a geometric weight. When two bonds form an angle $\theta$, the number of possible bond configurations in 3D space is proportional to $\sin(\theta)$.

A common example of bending energy is the harmonic potential:

$U_b(\theta)=-\frac{1}{2}k \theta^2$

where $k$ is the stiffness constant in units of $k_BT$. For this potential, $\langle\cos(\theta)\rangle$ is given by:

$\frac{e^{\frac{3}{4 k}} \left(-2 \text{erf}\left(\frac{1}{\sqrt{k}}\right)+\text{erf}\left(\frac{1+i \pi  k}{\sqrt{k}}\right)+\text{erf}\left(\frac{1-i \pi  k}{\sqrt{k}}\right)\right)}{2 \left(-2 \text{erf}\left(\frac{1}{2 \sqrt{k}}\right)+\text{erf}\left(\frac{1+2 i \pi  k}{2 \sqrt{k}}\right)+\text{erf}\left(\frac{1-2 i \pi  k}{2 \sqrt{k}}\right)\right)}$

For large values of $k$ ($k\ge 6$), $\langle\cos(\theta)\rangle$ approaches the Langevin function, $L(k)$.

A Gaussian chain, when averaged over all conformations and orientations, exhibits isotropic behavior. This means it can be treated as a sphere with a radius equal to its radius of gyration, $R_g$. However, an exercise in the second chapter of Rubinstein’s Polymer Physics demonstrates that $R_g^2$ becomes asymmetric when the coordinate frame is aligned with the end-to-end vector, $\mathbf{R}_{ee}$. Specifically, if $\mathbf{R}_{ee}$ is set as the x-axis, i.e., $\mathbf{R}_{ee}=(R,0,0)^T$, the three components of $R_g^2$ are $\frac{Nb^2}{36}$, $\frac{Nb^2}{36}$ in the y and z directions, and $\frac{Nb^2}{36}+\frac{R^2}{12}$ in the x direction. Here’s a simple proof:

Consider a Gaussian chain fixed between $(0,0,0)^T$ and $\mathbf{R}_{ee}$. This forms a Brownian bridge, with a multivariate Gaussian distribution centered at $\frac{i}{N}\mathbf{R}$ and a variance of $\frac{i(N-i)}{N}b^2$. The proof is straightforward:

$P_{0\mathbf{R}}(\mathbf{r},n)=\frac{G(\mathbf{r},0,n)G(\mathbf{R},\mathbf{r},N-n)}{G(0,\mathbf{R},N)}$

where $G(a,b,n)$ represents the distribution of a Gaussian chain with segment length $n$ and ends at points $a$ and $b$. This expression represents the probability of a length-$n$ Gaussian chain starting at 0 and ending at $\mathbf{r}$, connected to another length-$(N-n)$ Gaussian chain starting at $\mathbf{R}$ and ending at $\mathbf{r}$. The entire chain is a Gaussian chain with length $N$ and end-to-end vector $\mathbf{R}_{ee}=\mathbf{R}$. It’s easy to show that the distribution of this chain:

$P_{0\mathbf{R}}(\mathbf{r},n)=G\left(\mathbf{r}-\frac{n}{N}\mathbf{R}, 0,\frac{n(N-n)}{N}\right)$

is equivalent to a Gaussian chain segment with ends at $\mathbf{r}$ and $\frac{n}{N}\mathbf{R}$, and an equivalent length of $\frac{n(N-n)}{N}$. The $R_g^2$ is then:

$\begin{align}R_g^2=&\frac{1}{2N^2}\int_0^N\langle(\mathbf{r}_i-\mathbf{r_j})^2\rangle\mathrm{d}i\mathrm{d}j\\
=&\frac{1}{N^2}\int_0^N\int_j^N \frac{(i-j)^2 R^2}{N^2}+\frac{(i-j)(N-(i-j))}{N}b^2\mathrm{d}i\mathrm{d}j\\ =&\frac{R^2+Nb^2}{12}\end{align}$

In this equivalent method, $\langle (\mathbf{r}_i-\mathbf{r}_j)^2\rangle$ is interpreted as $(\frac{i}{N}-\frac{j}{N})^2R^2+\frac{|i-j|(N-|i-j|)}{N}$ from the equivalent Gaussian chain. This is because, despite the chain being ‘fixed,’ it remains a Gaussian chain, implying translation invariance. Therefore, $\langle (\mathbf{r}_i-\mathbf{r}_j)^2\rangle$ depends only on $|i-j|$, and $P_{0\mathbf{R}}(\mathbf{r},n)$ provides the probability of segment $\mathbf{r}_n-\mathbf{r}_0$, which represents any $n$-segment on the Gaussian chain. Directly calculating $P(\mathbf{r}_i-\mathbf{r}_j)$ from the convolution of $P_{0\mathbf{R}}$ is incorrect because $\mathbf{r}_i$ and $\mathbf{r}_j$ are correlated. Convolution of $P_{0\mathbf{R}}$ introduces dependence on $i$ and $j$ in $\langle (\mathbf{r}_i-\mathbf{r}_j)^2\rangle$, violating the translation invariance of the chain.

If $R^2=Nb^2$, we find that $R_g^2=\frac{1}{6}Nb^2$. Furthermore, we have $\frac{Nb^2}{36}+\frac{R^2_{x,y,z}}{12}=\frac{Nb^2+R^2}{36}$ for each dimension. Specifically, if $\mathbf{R}_{ee}=(R,0,0)^T$, we obtain $\frac{Nb^2}{36}$ in the y and z directions and $\frac{Nb^2}{36}+\frac{R^2}{12}$ in the x-direction.

Here’s a simple simulation code:

rg2_ree = rg2 = 0
for _ in  range(5000):  
    ch = np.random.normal(size=(1000,3))  
    ch = np.append(np.zeros((1,3)),np.cumsum(ch, axis=0),axis=0)  
    ch = ch - ch.mean(axis=0)  
    ree = ch[-1]-ch[0]  
    ree = ree/np.linalg.norm(ree)  
    rg2_ree += np.sum(ch.dot(ree)**2)  
    rg2 += np.sum(ch**2)
rg2_ree/rg2
# Result is 0.666...

1. This exercise focuses on the Flory characteristic ratio from Rubinstein’s textbook on Polymer Physics.
2. Definition of Flory Characteristic Ratio:
   
   
   $C_\infty:=\lim_{n\to\infty}\frac{\langle \mathbf{R}^2 \rangle}{nb^2}$
   
   
   where $\langle \mathbf{R}^2 \rangle$ represents the mean square end-to-end vector, $b$ is the Kuhn length of the polymer, and $n$ is the corresponding chain length.
   
3. Calculation of $\langle \mathbf{R}^2\rangle$:
   
   
   Let $\mathbf{r}_i$ denote the bond vector between the $i$th and $(i-1)$th particle, and let $\mathbf{r}_1=\mathbf{R}_1-\mathbf{R}_0\equiv0$. For an $n$-bead chain:
   
   
   $$\begin{align}\langle \mathbf{R}^2 \rangle&=\left\langle\left(\sum_{i=1}^n \mathbf{r}_i^2\right)\cdot\left(\sum_{j=1}^n \mathbf{r}_j^2\right)\right\rangle\\
   
   &=\sum_i\sum_j\langle \mathbf{r}_i\cdot\mathbf{r}_j\rangle\end{align}$$
   
4. Calculation of $\langle\mathbf{r}_i\cdot\mathbf{r}_j\rangle$:
   
   
   Consider a local coordinate frame where $\mathbf{r}_i=(b,0,0)^T$. Then, $\mathbf{r}_{i+1}$ can be viewed as a rotation of $\mathbf{r}_i$ by $\theta$ about the $z$-axis (bond angle) and by $\phi_{i}$ about the $x$-axis (torsion). The rotation matrix is:
   
   
   $$\mathbf{T}_i:=\left(
   
   \begin{array}{ccc}
   
   \cos (\theta ) & -\sin (\theta ) & 0 \\
   
   \sin (\theta ) \cos (\phi_i) & \cos (\theta ) \cos (\phi_i) & -\sin (\phi_i) \\
   
   \sin (\theta ) \sin (\phi_i ) & \cos (\theta ) \sin (\phi_i ) & \cos (\phi_i ) \\
   
   \end{array}\right)$$
   
   
   or $\pi-\theta$, $\pi – \phi$ based on the definitions of bond angle $\theta$ and torsion angle $\phi$. In each reference frame $R_i$, the bond vector is $(b,0,0)^T$, and in $R_i$, $\mathbf{r}_{i+1}=\mathbf{T}_i\cdot\mathbf{r}_i=(b \cos (\theta ),b \sin (\theta ) \cos (\phi_i ),b \sin (\theta ) \sin (\phi_i ))^T$.
   
   Now, consider the simple case of $\mathbf{r}_i\cdot\mathbf{r}_{i+2}$. $\mathbf{r}_{i+2}$ is calculated in the $(i+1)$th reference frame $R_{i+1}$ where $\mathbf{r}_{i+1}$ is $(b,0,0)^T$ and $\mathbf{r}_{i+2}=\mathbf{T}_{i+1}\cdot\mathbf{r}_{i+1}$. To perform the inner product, we need to rotate one of the vectors into the other frame. This involves rotating $\mathbf{r}_i=(b,0,0)^T$ from $R_i$ to $R_{i+1}$ or vice versa.
   
   Note that $\mathbf{r}_{i+1}$ in the $R_i$ frame is $(b \cos (\theta ),b \sin (\theta ) \cos (\phi_i ),b \sin (\theta ) \sin (\phi_i ))^T$ and $(b,0,0)^T$ in the $R_{i+1}$ frame. Therefore, to transform a vector from $R_i$ to $R_{i+1}$, we need a transformation matrix that converts $(b \cos (\theta ),b \sin (\theta ) \cos (\phi_i ),b \sin (\theta ) \sin (\phi_i ))^T$ to $(b,0,0)^T$. This matrix is $\mathbf{T}_i^{-1}$, or $\mathbf{T}_i^T$, because in frame $R_i$, we obtain $\mathbf{r}_{i+1}$ using $\mathbf{T}_i\cdot(b,0,0)^T$, and $\mathbf{T}_i$ is an orthonormal rotation matrix.
   
   By transforming $\mathbf{r}_i$ into $R_{i+1}$, the inner product $\mathbf{r}_i\cdot\mathbf{r}_{i+2}$ becomes:
   
   
   $$\langle\mathbf{T}_i^T\mathbf{r}_i,\mathbf{T}_{i+1}\mathbf{r}_{i+1}\rangle=$$
   $$\langle\mathbf{r}_i,\mathbf{T}_i\mathbf{T}_{i+1}\mathbf{r}_{i+1}\rangle=(b,0,0)\mathbf{T}_i\mathbf{T}_{i+1}(b,0,0)^T$$
   
   By induction, $\mathbf{r}_i\cdot\mathbf{r}_j=(b,0,0)\mathbf{T}_i\mathbf{T}_{i+1}\mathbf{T}_{i+2}\cdots\mathbf{T}_{j-1}(b,0,0)^T$. Since $\theta$ is the same for all bonds and $\phi_{1,2,\dots,n}$ are independently distributed, let $\mathbf{T}$ represent the average transformation matrix of ${\mathbf{T}_i}$:
   
   
   $$\begin{align}\mathbf{T}&:=\frac{\int_{-\pi}^{\pi}\mathbf{T}_i(\phi)\exp{(-U(\phi)/k_BT)}\mathrm{d}\phi}{\int_{-\pi}^{\pi}\exp{(-U(\phi)/k_BT)}\mathrm{d}\phi}\\
   
   &=\left(
   
   \begin{array}{ccc}
   
   \cos (\theta ) & -\sin (\theta ) & 0 \\
   
   \sin (\theta ) \langle\cos (\phi) \rangle & \cos (\theta )\langle \cos (\phi)\rangle &0 \\
   
   0& 0& \langle\cos (\phi)\rangle \\
   
   \end{array}\right)\end{align}$$
   
   Therefore, we have $\langle\mathbf{r}_i\cdot\mathbf{r}_j\rangle=(b,0,0)\mathbf{T}^{|j-i|}(b,0,0^T)=b^2(\mathbf{T}^{|j-i|})_{11}$. Every $\sin(\phi)$ term vanishes in the integral because $U$ is even.
   
5. Flory Characteristic Ratio:
   
   
   $$\begin{align}C_\infty&=\lim_{n\to\infty}\frac{1}{nb^2}\sum_i\sum_j \mathbf{r}_i\cdot\mathbf{r}_j\\
   
   &=\lim_{n\to\infty}\frac{1}{n}(\sum_i\sum_j T^{|j-i|})_{11}\\
   
   &=\lim_{n\to\infty}\left(-\frac{1}{n}\frac{-2 \mathbf{T}^{n+1}+\mathbf{T}^2 n+2 \mathbf{T}-n\mathbf{I}}{(\mathbf{I}-\mathbf{T})^2}\right)_{11}\\
   
   &=\left(\frac{\mathbf{I}+\mathbf{T}}{\mathbf{I}-\mathbf{T}}\right)_{11}\\
   
   &=\frac{1+\cos(\theta)}{1-\cos(\theta)}\frac{1-\langle\cos(\phi)\rangle}{1+\langle\cos(\phi)\rangle}\end{align}$$
   
   Here, $\frac{\bullet}{n}=0$ and $\mathbf{T}^n=0$ as $n\to\infty$ because $\mathbf{T}$ is constant, and the correlation between two beads approaches zero as the distance between them increases.