This Blog is Untitled

A circulant matrix is defined as:

$C=\begin{bmatrix}c_{0}&c_{n-1}&\dots &c_{2}&c_{1}\\c_{1}&c_{0}&c_{n-1}&&c_{2}\\\vdots &c_{1}&c_{0}&\ddots &\vdots \\c_{n-2}&&\ddots &\ddots &c_{n-1}\\c_{n-1}&c_{n-2}&\dots &c_{1}&c_{0}\\\end{bmatrix}$

where $C_{j, k}=c_{j-k \mod n}$. The $k$-th eigenvalue $\lambda_k$ and eigenvector $x_k$ satisfy $C\cdot x_k=\lambda_k x_k$, which can be expressed as $n$ equations:

$\sum_{j=0}^{m-1}c_{m-j}x_j+\sum_{j=m}^{n-1}c_{n-j+m}x_j=\lambda_k x_m\quad m=0,1,\dots,n-1$

with $c_n=c_0$, where $x_m$ is the $m$-th component of the eigenvector $x_k$. By changing the dummy summing variables ($j\to m-j$ and $j\to n-j+m$), we obtain:

$\sum_{j=1}^{m}c_j x_{m-j} +\sum_{j=m+1}^{n}c_j x_{n+m-j}=\lambda_k x_m$

with $m=0,1,2,\dots,n-1$. We can “guess” a solution where $x_j=\omega^j$, which transforms the equation into:

$\begin{align}&\sum_{j=1}^{m}c_j \omega^{m-j} +\sum_{j=m+1}^{n}c_j \omega^{n+m-j}=\lambda_k \omega^m\\ \leftrightarrow &\sum_{j=1}^{m}c_j \omega^{-j} +\omega^{n}\sum_{j=m+1}^{n}c_j \omega^{-j}=\lambda_k\end{align}$

Let $\omega$ be one of the $n$-th roots of unity, i.e., $\omega = \exp(-2\pi\mathbb{i}/n)$, then $\omega^{n}=1$. This leads to the eigenvalue:

$\lambda = \sum_{j=0}^{n-1}\omega^{-j} c_j$

with the corresponding eigenvector:

$x= (1, \exp(2\pi\mathbb{i}/n), \exp(2\pi\mathbb{i}/n)^2,\dots,\exp(2\pi\mathbb{i}/n)^{n-1})^T$

The $k$-th eigenvalue is generated from $\omega^{-k} = \exp(2\pi k\mathbb{i}/n)$ with $k\in [0,n-1]$, resulting in the $k$-th eigenvalue and eigenvector:

$\lambda_k = \sum_{j=0}^{n-1}c_j \exp(2\pi k\mathbb{i}/n)$

and

$x_k = (1, \exp(2\pi k\mathbb{i}/n), \exp(2\pi k\mathbb{i}/n)^2,\dots,\exp(2\pi k\mathbb{i}/n)^{n-1})^T$

The eigenspace is simply the DFT matrix, and all circulant matrices share the same eigenspace. It is easy to verify that circulant matrices possess the following properties:

If $A$ and $B$ are circulant matrices, then:

1. $AB=BA=W^\ast\Gamma W$ where $W$ is the DFT matrix and $\Gamma=\Gamma_A\Gamma_B$, with $\Gamma_i$ representing the diagonal matrix consisting of the eigenvalues of $i$; $\Gamma_A = WAW^\ast$.
2. $B+A=A+B=W^\ast\Omega W$, where $\Omega=\Gamma_A+\Gamma_B$.
3. If $\mathrm{det}(A)\ne 0$, then $A^{-1}=W^\ast \Gamma_A^{-1}W$.

The proof is straightforward:

1. $AB=W^\ast \Gamma_AW W^\ast\Gamma_BW=$
   $W^\ast\Gamma_A\Gamma_BW=W^\ast \Gamma_B\Gamma_AW=BA$
2. $W(A+B)W^\ast=\Gamma_A+\Gamma_B$.
3. $AW^\ast \Gamma_A^{-1}W=W^\ast \Gamma_AW W^\ast \Gamma_A^{-1}W=\mathbb{I}$

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