# Bethe Approximation for a Semi-flexible Polymer Chain

## Abstract

We present a Bethe approximation to study lattice models of linear polymers. The approach is variational in nature and based on the cluster variation method (CVM). We focus on a model with $`(i)`$ a nearest neighbor attractive energy $`ϵ_v`$ between pair of non–bonded monomers, $`(ii)`$ a bending energy $`ϵ_h`$ for each pair of successive chain segments which are not collinear. We determine the phase diagram of the system as a function of the reduced temperature $`t=\frac{T}{ϵ_v}`$ and of the parameter $`x=\frac{ϵ_h}{ϵ_v}`$. We find two different qualitative behaviors, on varying $`t`$. For small values of $`x`$ the system undergoes a $`\theta `$ collapse from an extended coil to a compact globule; subsequently, on decreasing further $`t`$, there is a first order transition to an anisotropic phase, characterized by global orientational order. For sufficiently large values of $`x`$, instead, there is directly a first order transition from the coil to the orientational ordered phase. Our results are in good agreement with previous Monte Carlo simulations and contradict in some aspects mean–field theory. In the limit of Hamiltonian walks, our approximation recovers results of the Flory-Huggins theory for polymer melting. PACS numbers: 05.70.Fh, 36.20.Ey, 64.60.Cn   

The configurational statistics of a long, linear polymer in solution has often been modeled by a self-avoiding walk (SAW) on a lattice . The self–avoiding constraint takes into account excluded–volume effects. Attractive Van der Waals interactions between monomers are also generally considered. They are included by assigning a negative energy $`ϵ_v`$ to each pair of nearest–neighbor monomers on the lattice, provided they are not consecutive along the chain. These interactions become relevant at low enough temperature, causing a collapse transition of the polymer. The transition point is called $`\theta `$ point and it separates a high temperature expanded structure from a low temperature compact globule.

Polymers with a local stiffness have also been considered, by including a bending energy $`ϵ_h`$ which favors straight segments of the chain. The semi-flexible model has attracted much interest in the low temperature phase and in particular in the limit of Hamiltonian walks, where the path is forced to visit all sites of the lattice. In this case, it is believed to describe the melting of polymers chains (see also ). The system undergoes a phase transition between a disordered (liquid) phase and an ordered (solid) phase, the latter being characterized by anisotropic orientational order. More recently, the semi-flexible model has attracted renewed interest because of the possible connection with the protein folding problem. In this spirit, each link of the walk represents an $`\alpha `$–helical turn (ca. $`3`$ amino acids) and the curvature term mimics the tendency to form secondary structures. The attractive energy between monomers models the hydrophobic effect, which is supposed to be the main driving force for the folding transition.

In the present work we consider a Bethe approximation for lattice homopolymers. Our aim is twofold: first we introduce a new method to deal with SAW problems, second we present results concerning the phase diagram of an isolated semi-flexible chain. In the approach to the problem we have followed the cluster variation method (CVM). This is a closed form approximation, which is known to give excellent results for the phase diagram of spin systems. The approximation scheme is determined by the largest clusters of sites which are treated exactly. The CVM allows us to write an approximate expression for the free energy of the system, as a function of the probability of occurrence of all possible configurations of the basic cluster. This free energy has then to be minimized, subject to consistency conditions on the distribution variables. The pair approximation considers nearest-neighbor pair of lattice sites and it corresponds to the Bethe approximation.

We represent a polymer as a SAW on a $`d`$–dimensional hyper-cubic lattice with $`V=L^d`$ sites. Pairs of non–consecutive vertices along the chain interact through an attractive nearest–neighbor potential $`ϵ_v`$. Stiffness of the chain is incorporated by attributing an energy penalty $`ϵ_h`$ to each turn (corner) of the walk. Let $`T`$ be the absolute temperature and $`\beta =\frac{1}{\kappa _BT}`$ the inverse temperature. In the following we will adopt the notation $`\omega =\beta ϵ_v`$, $`t=\frac{T}{ϵ_v}`$ and $`x=\frac{ϵ_h}{ϵ_v}`$. The partition function of the system is

$$Z_N=\underset{\{SAW\}}{}e^{\omega (N_{con}(𝒮)xN_{cor}(𝒮))}$$
(1)   
where $`\{SAW\}`$ denotes the ensemble of all $`N`$–step SAW; $`N_{con}(𝒮)`$ and $`N_{cor}(𝒮)`$ are respectively the number of contacts and corners in walk $`𝒮`$. Introducing a monomer fugacity $`z`$, the grand canonical partition function reads therefore

$$𝒵=\underset{N=1}{\overset{\mathrm{}}{}}z^NZ_N$$
(2)   
where the sum is over all possible lengths $`N`$ of the walk.

The lattice model (1) (or equivalently (2)) has been the object of recent investigations. A mean–field analysis predicts a $`\theta `$–collapse transition at a temperature $`t_\theta `$ independent of $`x`$. Another transition should occur at lower temperature. It is a discontinuous melting transition from a disordered globule to an ordered “crystalline” phase. The melting temperature $`t_m`$ increases with $`x`$, although in mean-field approximation $`t_m<t_\theta `$ for any value of $`x`$. This picture has been partly contradicted by heuristic arguments and by Monte Carlo simulations. Indeed $`t_\theta `$ appears to slightly increase with $`x`$. Most importantly, $`t_m`$ seems to grow with stiffness and it does not reach a finite asymptotic value. This implies that the line of the melting transition hits the line of the $`\theta `$–collapse. For sufficiently high values of $`x`$, there is therefore a direct first order transition from the open coil to the ordered phase. The triple point is found to be approximately at $`x13`$.

The formulation of the CVM given by An  is a particularly convenient starting point for our analysis. We first illustrate the method by treating in some detail the simpler case of zero stiffness ($`x=0`$). Let $`s_i`$ and $`p_i`$ be the distribution variables assigned respectively to each site and pair configurations. One should in principle distinguish among all the possible configurations which are not related by symmetry operations. In $`d=2`$ there are, for instance, $`3`$ single site and $`11`$ pair independent configurations. In fact we have verified that in the particular case $`x=0`$ they can be grouped into a smaller number of non–equivalent classes. These are determined only by the following conditions: (a) the site is visited by the path, (b) the nearest–neighbor pair is joined directly by the path. The independent configurations are reported schematically in fig. 1, together with their multiplicity of occurrence. Following the notation of fig. 1, the free energy of the system (2), in the pair approximation, can be written

$`{\displaystyle \frac{\beta F}{V}}`$ $`=`$ $`{\displaystyle \frac{q(q1)}{2}}\mathrm{ln}zs_1{\displaystyle \frac{q(q1)^2(q2)^2}{8}}\omega p_2+`$ (4)   
$`(1q)\left({\displaystyle \underset{i=1}{\overset{2}{}}}m_s(i)s_i\mathrm{ln}s_i\right)+{\displaystyle \frac{q}{2}}\left({\displaystyle \underset{i=1}{\overset{4}{}}}m_p(i)p_i\mathrm{ln}p_i\right)`$   
where $`q=2d`$ is the coordination number of the lattice and $`m_s(i)`$ and $`m_p(i)`$ stand respectively for the multiplicity of site and pair configurations. Normalization of the distributions and consistency conditions on the probability variables require respectively

$$s_2=1\frac{q(q1)}{2}s_1$$
(5)   
and

$`p_1`$ $`=`$ $`{\displaystyle \frac{s_1}{(q1)}}`$ (6)   
$`p_2`$ $`=`$ $`{\displaystyle \frac{2s_1}{(q1)(q2)}}{\displaystyle \frac{2}{(q1)(q2)}}p_3`$ (7)   
$`p_4`$ $`=`$ $`1{\displaystyle \frac{q(q1)}{2}}s_1{\displaystyle \frac{(q1)(q2)}{2}}p_3`$ (8)   
This leaves us with only two variational parameters, e.g. $`s_1`$ and $`p_3`$. Substituting (5) and (7) into (4) and minimizing with respect to $`s_1`$ and $`p_3`$ we obtain the stable phase at a given $`z`$ and $`\omega `$. We report the complete phase diagram for $`d=3`$ in fig. 2. The polymer is a critical system along the transition line $`z_c(\omega )`$. This line separates a chain with zero density ($`s_1=0`$ for $`z<z_c(\omega )`$) from a chain with finite density ($`s_10`$ for $`z>z_c(\omega )`$). The continuous line represents a second order transition and the average number of monomers diverges with a power law as $`z`$ tends to $`z_c`$. The broken line is instead a first order transition and the density of monomers makes a finite jump at $`z_c`$. The cross denotes the tricritical point and it corresponds to the $`\theta `$ point. In the case of pure SAW ($`\omega =0`$) the connectivity constant is $`\mu =z_c^1=2d1`$. This result would have been expected by studying SAW on a Bethe lattice and it should be compared, for instance, with $`\mu 2.64`$ and $`\mu 4.68`$, obtained from exact enumerations respectively in $`d=2`$ and $`d=3`$. In our framework $`z_c`$ does not depend on $`\omega `$, as long as $`\omega <\omega _\theta `$. This is certainly an artifact of the approximation. It can be ascribed to the fact that, after minimization in $`p_3`$, $`p_2s_1^2`$ in the limit of $`s_1`$ going to zero. As a consequence, there is no term in (4) proportional to $`\omega `$ which contributes in locating the minimum of the free energy around $`s_1=0`$. Nonetheless the estimates we obtain for the $`\theta `$ point $`\omega _\theta ^{(B)}`$ are a better approximation to the available numerical values, with respect to mean field theory ($`\omega _\theta ^{(MF)}=\frac{1}{2d}`$): in $`d=2`$, $`\omega _\theta 0.665`$, $`\omega _\theta ^{(B)}0.4055`$ and $`\omega _\theta ^{(MF)}=0.25`$; in $`d=3`$, $`\omega _\theta 0.275`$, $`\omega _\theta ^{(B)}0.2231`$ and $`\omega _\theta ^{(MF)}=0.1667`$.

In the general case $`x0`$ equations (4), (5) and (7) must be generalized to include the curvature energy and the possibility of an anisotropic phase. In this case it is not possible to group configurations as in fig. 1 and one has to face a complex minimization problem. An efficient way of doing it numerically is by mean of the natural iteration method. The resulting phase diagram is reported in fig.3, as a function of $`x`$ and $`t`$ for $`d=3`$. The fugacity $`z`$ is fixed to its critical value $`z_c(x,t)`$. This condition assures we are studying a polymer in the limit of infinite chain length ($`N\mathrm{}`$). We find three different phases: an open coil, a compact globule and an ordered crystal. In our approximation, the latter is just the ground state of the polymer, having all links perfectly aligned. This is known to be not completely correct, as it has been shown rigorously that, for instance, in the case of Hamiltonian walks the entropy strictly vanishes only in the limit $`T0`$. The $`\theta `$–collapse line between the coil and the globule appears to be independent of $`x`$. On the other hand the discontinuous melting transition tends to infinity with $`x`$. Beyond the triple point at $`x8.8`$, there is directly a first order transition from the coil to the solid.

The limit $`T0`$ (or $`z\mathrm{}`$) corresponds to Hamiltonian walks. In this limit, walks are space filling and configurations with vacancies do not contribute to the partition function (2). Also, the attractive nearest–neighbor potential $`ϵ_v`$ plays no role as there are precisely $`(d1)`$ contacts per monomer. The resulting model is the so called Flory model of polymer melting. In this case, we obtain a minimum for the free energy in the disordered, compact phase which coincides with the free energy estimated by using the Flory-Huggins approximation. The latter was originally derived from combinatorial arguments. Neglecting the constant contribution proportional to $`ϵ_v`$, the analytical expression reads

$$\frac{\beta F_{FH}}{V}=\mathrm{ln}\left[\frac{\left(1\frac{2}{q}\right)^{\left(\frac{q}{2}1\right)}}{1+(q2)\mathrm{exp}\left(\omega x\right)}\right]$$
(9)   
At low temperature $`F_{FH}`$ competes for stability with the local minimum associated to the ordered phase, which has strictly $`F_O=0`$ in our approach. A first order phase transition takes therefore place at

$$t_m=\frac{x}{\mathrm{ln}\left[\frac{q2}{\left(1\frac{2}{q}\right)^{\left(\frac{q}{2}1\right)}1}\right]}$$
(10)   
In particular for $`d=3`$ ($`q=6`$) we have $`t_m/x=\left(\mathrm{ln}\left(\frac{16}{5}\right)\right)^10.86`$, which corresponds to the slope of the globule-solid transition line of fig. 3, in the limit $`x0`$. This value is slightly larger than the analogous mean-field estimate, $`t_m/x0.58`$.

In the case $`x=0`$ each HW is equally weighted and the total number of paths is believed to scale as $`𝒩_{HW}\mu _H^N`$. From (9) we therefore obtain

$$\mu _H^{(B)}=q\frac{q1}{q2}(1\frac{2}{q})^{\frac{q}{2}}$$
(11)   
A mean-field approach, in very good agreement with numerical estimates, yields

$$\mu _H^{(MF)}=\frac{q}{e}$$
(12)   
It is interesting to note that expression (11) correctly predicts $`\mu _H=1`$ for $`d=1`$, contrary to (12). In order to evaluate corrections to mean–field theory in powers of $`1/q`$, we have formulated the Hamiltonian walk problem through spin variables. This is in close analogy to the De Gennes theorem for SAW. By using a suitable high–temperature expansion we have then been able to calculate the coefficients of the series up to the third order, extending of one order previous results. We find

$$\mu _H=\frac{q}{e}(1+\frac{1}{6q^2}\frac{2}{3q^3}+\mathrm{})$$
(13)   
In the limit $`q\mathrm{}`$ expression (11) specializes to

$$\mu _H^{(B)}=\frac{q}{e}(1+\frac{1}{6q^2}+\frac{1}{3q^3}+\mathrm{})$$
(14)   
The Bethe approximation therefore correctly recovers results of the expansion (13) up to order $`O(\frac{1}{q^2})`$.

In summary, we have presented a Bethe approximation for lattice models of linear homopolymers. The method constitutes a substantial improvement with respect to mean–field theory. Indeed it produces a phase diagram for a semi-flexible polymer chain which is in good agreement with Monte Carlo simulations. In particular, we find a triple point where the $`\theta `$ collapse line and the melting transition line meets. In the limit of Hamiltonian walks it recovers results of the Flory-Huggins theory for polymer melting, whose variational nature appears in a transparent shape within our framework. It has the advantage of not requiring any spin or field theoretical representation, rather it relies directly on the configurations of the system. This last consideration suggests the scheme is more general and suitable to be applied to other geometrical problems, as, for instance, branched polymers and self-avoiding surfaces (see  and references therein). It is plausible that the accuracy of the method can be systematically refined according to the cluster variation method, in analogy with spin systems. We expect some of the inaccurate features of the pair approximation could be removed by considering larger basic clusters.

<sup>a</sup> Electronic address: lise@sissa.it   
<sup>b</sup> Electronic address: maritan@sissa.it   
<sup>c</sup> Electronic address: alex@athena.polito.it