Multidimensional integral to compute the end-to-end distribution of a FENE ideal chain

Question

How do I evaluate the mean end-to-end squared distance of a FENE ideal chain at fixed inverse temperature $\beta$ in the canonical ensemble?

This quantity is defined as the mean value $$\left<\left(\sum_{i=1}^{N-1} \vec r_i\right)\cdot\left(\sum_{j=1}^{N-1} \vec r_j\right)\right>$$

Under the canonical distribution $$ \exp\left(-\beta \sum_{i=1}^NH(\vec r_i)\right)\text{ where } H(\vec r_i) = - \ln\left(1 - \left(\frac{|\vec r_i| - r_0}{\Delta}\right)^2\right)\,. $$

Therefore, it would be $$ \frac{ \int \prod_{i=1}^{N-1}d\vec r_i\left(1 - \left(\frac{|\vec r_i| - r_0}{\Delta}\right)^2\right)^\beta \left(\sum_{i=1}^{N-1} \vec r_i\right)\cdot\left(\sum_{j=1}^{N-1} \vec r_j\right)} {\int \prod_{i=1}^{N-1}d\vec r_i\left(1 - \left(\frac{|\vec r_i| - r_0}{\Delta}\right)^2\right)^\beta} $$ where the integral is to be performed over all the values of $\vec r_i$ that give a positive argument of the logarithm. Edit: I realised that the integral that I wanted to solve was not the one written below, but the one written here above. The following is the question in its first (wrong) version, which is the one to which Yarchik responded.

How do I evaluate the quantity $\left.\frac{d}{d\alpha} \ln Z_\alpha\right|_{\alpha=0}$ for N=30 in Mathematica, when

$$ Z_\alpha = V\int e^{-\beta [H(\vec r_1)+H(\vec r_2-\vec r_1)+...+H(\vec r_{N-1}-\vec r_{N-2})]+\alpha|r_{N-1}|^2}d\vec r_1d\vec r_2\cdots d\vec r_{N-1} $$ and $$ H(\vec r) =- \ln\left(1-16(|\vec r|-1)^2\right) ? $$ The integration is extended to all the domain in which the integral is defined (that is, the argument of the function $H$ must be a vector with length in [3/4,5/4]).

I know how to declare the function H, but how can I nest 30 integrate functions without writing them explicitly? I don't want to rewrite my code from scratch if I change the value of N I wish to consider. So far I have the following (H is expressed in spherical coordinates, so that r is just the distance between the particles):

a = 1; d = 0.25; r0 = 1
H[r, r0, d] = -a*Log[1 - (r - r0)^2/d^2]
U = FullSimplify[-D[Log[Z[r0, d, b]], b]]

but I do not know how I can write Z.

For the record, I need this to compute the end to end distance distribution of an ideal chain of 30 atoms linked by FENE springs. $H$ is the hamiltonian associated with each of the 29 springs, and Z is the partition function of my system. The derivative I want to compute is a way of computing the mean of the parenthesis that multiplies $\alpha$.

Edit: thanks to Yarchik for noting that the domain of H was wrong. It should be right now, but in case there's an error I want to integrate over all the vectors for which H is defined.

Edit: I also realised that I was writing the end to end distance in a wrong way. As Yarchik pointed out, the way I was writing it I was computing the position of the center of mass, while now I'm actually computing the average end to end squared distance.

Answer 1

I believe it is not possible to solve this problem with Mathematica in the presented form. Partly because there inconsistencies in the formulation:

Plot[-Log[1 - 16 (Abs[x] - 1)^2], {x, 1 - 1/4, 1 + 1/4}]

We see that the potential is, in fact, defined on a larger domain [3/4,5/4]

An approximate solution can be numerically obtained using a Monte Carlo integration. Notice, the expectation value of the center-of-mass position has to be computed directly, rather then through the partition function. That is to answer the question on how to formulate it for Mathematica in $N$-invariant way. Below I would like to consider a special analytically solvable case as is also suggested by others. We take $N=2$, i.e. a chain consisting of a single particle at infinite temperature, i.e. $\beta=0$. In this case we deal with a 3d integral that in view of the spherical symmetry can be reduced to 1d.

  Z[a_] := 4 \[Pi] Integrate[Exp[a r^2] r^2, {r, 1-1/4, 1 + 1/4}] 

The answer is then

 Limit[D[Log[Z[a]], a], a -> 0]
 (*   4323/3920  *)

Analogically, the case of finite temperatures can be considered and analytically integrated for integer values of $\beta$:

 r = Table[{b, 
 Limit[
  D[Log[
  Integrate[
   Exp[a r^2] r^2 (1 - 16 (r - 1)^2)^b, {r, 1 - 1/4, 1 + 1/4}]], a], 
   a-> 0]}, {b, 0, 10}]
 fr = Interpolation[r]
 Plot[fr[b], {b, 0, 10}]

Finally, consider another extreme $N\rightarrow\infty$. Just for physical clarity it makes sense to redefine alpha such that the whole Hamiltonian has a form of interacting points with an extra potential acting on the center-of-mass of the chain. Notice, that in the case each particle is subject to this (quadratic) potential individually the problem is solvable.

Edit

Concerning the modified question on how to determine the distance to the end point $d$. An exact solution is probably not feasible, but a very good approximation can be suggested. First, modify the inter-particle potential $V(x)=-\ln(1-(x-1)^2/b^2)$. In the present setup $b=1/4$. What happens when $b\rightarrow 0$? The system is essentially a chain with fixed link-length ($=1$). In this case the end-to end distance is exactly $\sqrt{N}$. This follows from the fact that the average value $\langle \cos \theta\rangle=0$, where $\cos\theta$ is the angle between two neighbouring links. In the case when $1>b>0$ a very accurate result can be obtained using calculations in hands. Second figure shows the average link-length squared $\langle r^2\rangle$ as a function of the inverse temperature $\beta$. The end-to-end distance is therefore:

$d(\beta)\approx\sqrt{N\langle r^2\rangle(\beta)}.$

When temperature is zero, i.e., $\beta=\infty$ another exact result reads $d(\infty)=\sqrt{N}.$

Answer 2

First, let's rewrite the mean value as the sum of the mean values of the scalar products of the single particle relative positions:

\begin{align*}r^2_{e2e} & = \left<\left(\sum_{i=1}^{N-1} \vec r_i\right)\cdot\left(\sum_{j=1}^{N-1} \vec r_j\right)\right>\\ &=\sum_{i=1}^{N-1} \sum_{j=1}^{N-1} \left< \vec r_i\cdot \vec r_j\right>\\ & = (N-1) \left<|\vec r_i|^2\right>+\sum_{i<j}\left<\vec r_i\cdot \vec r_j\right>\\ & = (N-1) \left<|\vec r_i|^2\right>+\sum_{i<j}\left<\vec r_i\right>\cdot \left<\vec r_j\right>\\\end{align*} Where the last passage can be done since the different positions are uncoupled, i.e. evolve indipendently. The interesting thing is that both the sincgle vector means in the last passage are both zeroes, since each component yelds an odd function integrated on a symmetric domain. Therefore, we can focus on computing the first average in the last passage. Averaging over all the degrees of freedom different from $i$ yields a factor of one, therefore now the problem is simply to find the average distance between two particles connected by the FENE potential (that is, by using a single coordinate).

As is customary in statiscal mechanics, let us define the one oscillator partition sum (that is,the normalisation constant of the system with only two particles), and solve it in polar coordinates. Since $H$ is radial, we get: $$ Z = \int d\vec r_1\left(1 - \left(\frac{|\vec r_1| - r_0}{\Delta}\right)^2\right)^\beta = 4\pi\int^{r_0+\Delta}_{r_0-\Delta} dr \left(1 - \left(\frac{r - r_0}{\Delta}\right)^2\right)^\beta r^2\,. $$

Similarly we could write the numerator of the quantity that we wish to calculate, finally yielding:

$$ r^2_{e2e} = (N-1)\frac{4\pi\int^{r_0+\Delta}_{r_0-\Delta} dr \left(1 - \left(\frac{r - r_0}{\Delta}\right)^2\right)^\beta r^4}{4\pi\int^{r_0+\Delta}_{r_0-\Delta} dr \left(1 - \left(\frac{r - r_0}{\Delta}\right)^2\right)^\beta r^2} $$

This can be computed with the Mathematica function Nintegrate (for some reason, using Integrate instead will give problems for some values of $\beta$, like $\beta = 10$.

(* by this point, b, N, r0 and D must have already been defined*)
I[a_,b_]:=NIntegrate[(1 - (r - r0)^2/D^2)^b r^a, {r,r0-D,r0+D}]
( N - 1 ) I[4,b]/I[2,b]