Phase Transitions in Pressurised Semiflexible Polymer Rings
We propose and study a model for the equilibrium statistical mechanics of a pressurised semiflexible polymer ring. The Hamiltonian has a term which couples to the algebraic area of the ring and a term which accounts for bending (semiflexibility). The model allows for self-intersections. Using a combination of Monte Carlo simulations, Flory-type scaling theory, mean-field approximations and lattice enumeration techniques, we obtain a phase diagram in which collapsed and inflated phases are separated by a continuous transition. The scaling properties of the averaged area as a function of the number of units of the ring are derived. For large pressures, the asymptotic behaviour of the area is calculated for both discrete and lattice versions of the model. For small pressures, the area is obtained through a mapping onto the quantum mechanical problem of an electron moving in a magnetic field. The simulation data agree well with the analytic and mean-field results.
Fluid vesicles obtained via the self-assembly of amphiphilic molecules exhibit a variety of shapes in thermal equilibrium. Such shapes can be understood in terms of the energy minimising configurations of a curvature Hamiltonian, under the constraints of fixed enclosed volume and surface areaCanham (1970); Helfrich (1973); Evans (1974). Shape changes arise when solutions of the Euler-Lagrange equations representing distinct shapes exchange stability. However, the non-linearity of these equations, if no special symmetries are assumed, necessitates purely numerical approaches. Further, while the curvature modulus in bilayer lipid membrane systems is often large, so that thermal fluctuations about the minimum free energy structure may be ignored, the more general problem of understanding the thermodynamics of such shape transitions is a formidable oneSeifert (1997).
The two-dimensional version of the vesicle problem is a polymer ring of fixed contour length, whose enclosed area is constrained through a coupling to a pressure difference term . Leibler, Singh and Fisher (LSF) Leibler et al. (1987) performed a Monte Carlo and scaling study of two-dimensional vesicles, modelled as closed, planar, self-avoiding tethered chains, accounting for both pressure and bending rigidity. In this model, the ring polymer is obtained by connecting the centres of impenetrable particles of fixed radius with tethers of a fixed maximum length, while enforcing self-avoidance. LSF showed the existence of a phase transition at , separating a branched polymer phase for from an inflated phase for . At the transition point, the ring is described by a self avoiding polygon. Various fractal and non-fractal shapes that arise in these models have also been investigated Fisher (1989); Camacho and Fisher (1990).
Analytic studies of this class of models present many difficulties, arising principally from the self-avoidance constraint. Nevertheless, the relatively simple structure of the LSF model has stimulated a considerable body of work, largely in exact enumeration studies of lattice versions of the original continuum model and its variants van Rensburg (2000); Fisher et al. (1991); Bousquet-Mélou (1996); Guttmann and Lin (1992); Richard et al. (2001); Cardy (2001); Richard (2002); Richard et al. (2003); Rajesh and Dhar (2005). Most of these studies have concentrated on the behaviour of the system in the thermodynamic limit in the region . However, the case can exhibit interesting crossover behaviour for large but finite systems.
The consequences of relaxing the self-avoidance constraint were studied in Refs. Rudnick and Gaspari (1991); Gaspari et al. (1993); Marconi and Maritan (1993). In the models studied in these papers, the ring was allowed to intersect itself, with the pressure term coupled to the algebraic area Gaspari et al. (1993); Rudnick and Gaspari (1991) or to its square Marconi and Maritan (1993). The particles linked to form the polymer were coupled through harmonic springs Gaspari et al. (1993); Rudnick and Gaspari (1991), thus allowing for the extensibility of the chain. We shall refer to this model as the Extensible Self-Intersecting Ring (ESIR). The ESIR model can be solved exactly. The solution yields collapsed and inflated phases of the ring separated by a continuous phase transition that occurs at a critical value of an appropriately scaled pressureGaspari et al. (1993). However, the model has a major shortcoming in that the inflated phase is an unphysical one in which the ring expands to an infinite size. In a more realistic model, such an expansion would be limited by the finite size of individual link lengths.
The unphysical nature of the inflated phase in the ESIR model has been addressed in recent workHaleva and Diamant (2006), in which particles are joined by bonds of fixed length, as opposed to springs. The Hamiltonian has a term where the pressure couples to the algebraic area, as in the ESIR model. The transition survives as a continuous phase transition with mean-field exponents, separating collapsed and inflated regimes of the ring. We shall refer to this model as the Inextensible Self-Intersecting Ring (ISIR).
The model proposed in this paper incorporates a bending energy into the ISIR model along standard lines for semiflexible polymers. We retain the coupling of the signed pressure to the algebraic area noting, as argued in Haleva and Diamant (2006), that this difference, while vastly increasing the tractability of the problem, makes little difference to computations within the inflated phase.
The continuum problem we address is the following. If the polymer chain is specified by the curve , where is the arc-length along the curve, we consider the Hamiltonian
Here, is the infinitesimal (signed) area which must be integrated over the internal variable to obtain the total algebraic area. The quantity is the length of the polymer ring of units where is the size of the basic monomer. The continuum limit is taken such that and , keeping fixed. The parameter represents the pressure differential between the inside and the outside of the ring and is the bending rigidity of the chain. We measure energies in units of . The inextensibility constraint is imposed through
where is the unit tangent vector. When , this model reduces to the worm-like chain model, with configurations constrained by the closure requirement.
We will work with two discretized versions of the above Hamiltonian. The first has particles in the continuum connected through fixed length links, forming a ring polymer whose equilibrium configurations are constrained by the pressurisation and bending energy terms above, but where the ring can intersect itself at no energy cost. We will refer to this version of our model as the “discrete model”. The second is a square lattice version of the same problem with the particles constrained to lie on the vertices of the lattice. We will refer to this version as the “lattice model”. We discuss the differences and similarities between the two versions.
We use a combination of analytic and numerical methods to study these models: Flory type scaling theory for the scaling of the area as a function of pressure, Monte Carlo simulations for different pressures and bending rigidity, mean field approaches and exact enumeration.
In Fig. 1 we show typical configurations obtained from Monte Carlo simulations of the discrete model in four limits. These are configuration snapshots across the collapsed to inflated phase transition, for different values of the bending rigidity of the discrete model (), as the pressure is varied. Fig. 1(a) shows the collapsed phase for the case where the bending energy is zero, while Fig. 1(b) illustrates a typical ring configuration at an intermediate value of the bending rigidity, but still within the collapsed regime. In Fig. 1(c), we show a typical configuration close to the transition between collapsed and inflated phases. Last, Fig. 1(d) illustrates the fully inflated ring.
We summarise our main results below. We show that there is a continuous phase transition in the scaled pressure – bending rigidity () phase diagram, which separates a collapsed phase in which , from an inflated phase in which (see Fig. 2). The phase boundary for the discrete model is obtained as , where the ’s are modified Bessel functions. For the lattice model, the phase boundary is obtained as .
These results are obtained by first solving the case exactly and then incorporating the effects of a nonzero through a scaling argument. For the collapsed phase, the free energy for nonzero is calculated by the same method. In the inflated regime, we resort to mean field theories. We employ two types of mean-field theories: In the first, the inextensibility constraint is satisfied exactly but the closure condition is satisfied only on average. In the second, we impose the closure condition exactly but satisfy the inextensibility constraint only on average. The dependence of the area on for is calculated. The behaviour near the transition line is obtained through a Flory type scaling theory.
The rest of the paper is organised as follows. In Sec. II we define our models more precisely. Section III contains the details of the numerical methods used, including the Monte Carlo and exact enumeration algorithms. In Sec. IV, we discuss a Flory-type scaling theory valid for the semi-flexible case. Section V describes mean-field approaches to this problem: (a) a simple density-matrix based single-site mean-field approach which captures the properties of the inflated phase to very high accuracy but is inadequate for the collapsed phase and (b), a less accurate harmonic spring mean-field theory, which is capable of describing both collapsed and inflated phases. In Sec. VI, we discuss the behaviour around the critical point in greater detail. Our Sec. VII contains results for the asymptotic behaviour of the area as well as a description of the appropriate scaling function for the area in the lattice case, as a function of . Section VIII contains a summary and conclusions. In Appendix A we present a solution of the problem when by drawing an analogy with the quantum mechanical problem of the motion of an electron in a magnetic field.
Consider a closed chain of N monomers in two dimensions. Let the positions of the particle be denoted by the vector and the corresponding tangent vectors by For a closed ring, , or equivalently, . The algebraic or signed area area enclosed by the ring is given by
can be either positive or negative.
Coupling this algebraic area to pressure, we obtain the energy term,
The bending energy cost can then be written down following standard procedures as
where the bending rigidity of the discrete model is proportional to the continuum bending rigidity and is the unit vector in the direction of . The inextensibility condition is imposed through
Since the tangent vectors have unit norm, we can represent them as , where . In terms of these variables, the partition function is
On a square lattice, the model remains essentially the same except for restrictions on the angles . Now is allowed to take values such that all the particles are on the vertices of the square lattice.
Iii Numerical method
In this section, we describe the numerical methods used. For the discrete version of our model, we use Monte Carlo simulations (described in Sec. III.1) while for the lattice problem on the square lattice, we use an exact enumeration scheme (described in Sec. III.2). The analytic results we obtain for our model, described in later sections, provide useful benchmarks for the numerical work.
iii.1 Monte Carlo Simulations
The algorithm for the Monte Carlo simulation of the discrete model consists of two basic movesKoniaris (1994); van Rensurg et al. (1990): a single particle flip and a global flip. In the single particle flip, a particle is picked at random and reflected about the straight line joining its two neighbours (see Fig. 3(a)). The move is accepted using the standard Metropolis algorithm. Since the energy computation involves only nearby sites, the move is efficient and fast. In the global flip, two particles of the ring are chosen at random and the section of the ring between them is reflected about the line joining the two particles (see Fig. 3(b)). The energy calculation now involves particles and is thus computationally expensive. However, the global move is crucial to the study of the case where , since single particle moves alone are insufficient for equilibration in this case.
In the simulations, one Monte Carlo step is defined as one global move and single particle moves made by selecting at random particles to be updated. This step is then repeated until the system equilibrates. Thermodynamic quantities are measured from averages taken over independent configurations in equilibrium.
The initial configuration was chosen to be a regular N-sided polygon, but we verified that random configurations also gave the same results. We performed Monte Carlo simulations across a range of pressures for different values of J and system size. The system size varied from to . Typically each parameter value was run for Monte Carlo steps. We waited typically for steps for equilibration, averaging data over the remaining steps using independent configurations.
iii.2 Exact enumeration
We first describe the algorithm for the case . Consider a random walk starting from the origin and taking steps in one of the four possible directions. For each step in the positive (negative) -direction, we assign a weight (), where is the ordinate of the walker. Multiplying these weights, it is easy to check that the weight is for a closed walk enclosing an area .
Let be the weighted sum of all -step walks from to . It then obeys the recursion relation,
with the initial condition
Finally, gives the partition function of the ring polymer on a lattice.
For the semiflexible case, the recursion relation given above must be modified, since the ring is no longer a simple random walk but a walk with a one step memory. We convert it into a Markov process as follows. Let be the sum of weights of all walks reaching in steps but having been at at the previous step. These ’s are now a Markov process and depend only on ’s. The recursion relations are then straightforward to write down. Rather than give all the recursion relations, we provide a representative example
Similar recursion relations will hold for , and .
The partition function for the polymer problem can be expressed as a sum over areas and bends consistent with a given value of the area, i.e.,
where counts the number of closed paths of area in a walk of length N which have bends.
We count up to for different values of . The only limiting factor in going to larger values is computer memory.
Iv Flory-type Scaling Analysis
Flory type scaling theory provides a useful tool to capture the scaling behaviour of systems whose free energy reflects a competition between two or more terms. Such a scaling theory was proposed for the ISIR model in Ref. Haleva and Diamant (2006). A transition from a collapsed to an inflated state was predicted to occur at a critical value of the pressure, whose magnitude scaled with system size as . We show how these arguments may be extended to the semiflexible case, deriving expressions for the change in the critical point and scaling as a function of the bending rigidity.
The free energy consists of three terms describing (i) the entropy of the ring, (ii) the pressure differential and (iii) inextensibility of the bonds. When , these terms were argued to be , and for a ring of size Haleva and Diamant (2006), where for the second term it was assumed that the area scales as . With semiflexibility, we show that a similar scaling form holds except for dependent prefactors. Thus, the free energy takes the form
where we have defined , and and depend on .
It is easily seen that a system described by such a Flory theory undergoes a continuous transition when the sign of the term changes sign. This occurs at a critical scaled pressure which varies with as
When , then the area follows random walk statistics with . In this regime the term is not important. For length scales much larger than the persistence length, the problem is effectively one of a freely jointed ring, but with a suitably defined . Thus, we conclude that
where is a scaling function. The scaling function and can be determined by solving the case (see Appendix A). This gives
When , the scaling is determined by the term. Thus, . Thus,
To test this relation, we compare the Flory prediction with the enumeration results for the area in the lattice model. As can be seen from Fig. 5, there is good agreement for small values of but the data starts to deviate away from the predicted curve as increases.
When , the ring is in an inflated state, with the area . To obtain an accurate description of this regime, we would need to keep higher order terms such as and so on. One thus expects that the lattice and the discrete problems should differ considerably in this regime.
We now derive expressions for and in both the discrete and lattice cases. This is done by considering a semiflexible chain subjected to an external force. We obtain a perturbative solution for the partition function in the limit of small forces. From the partition function, we obtain the free energy of the ring. By comparing this with the form of Eq. (12), the values of and can be obtained.
iv.1 Discrete Case
Consider a semiflexible chain of N monomers. When the chain is pulled by a force , the partition function is given by
We work in the limit of small forces, treating the term exactly. We consider the term as a perturbation on the zeroth order partition function [ in Eq. (18)], given by
where is the modified Bessel function of the first kind of order . We then expand as a series in and average each term with respect to the zeroth order Hamiltonian. On computing the averages, the partition function is obtained as
where the coefficients and are given by
The ’s are modified Bessel functions of the first kind. Their dependence has been suppressed in the equation above.
iv.2 Lattice Case
For a lattice polygon, where each individual step can point only in four directions, we solve the problem of a semiflexible chain subject to an external force using the exact transfer matrix. The transfer matrix in this case is given by
We determine the largest eigenvalue up to order , and hence calculate the partition function:
We then follow the same procedure as for the discrete case, finding in terms of , inverting this equation to find , and finally using this expression to compute the free energy. We thus obtain
The expressions for and are then
V Mean Field Theory
In this section we present mean-field theories to calculate the dependence of area on pressure and bending rigidity. In Sec. V.1, we address the ISIR model (). The mean field theory presented in Haleva and Diamant (2006) performs poorly with respect to the Monte Carlo data when . Here, we present an improved variational mean field which reproduces the behaviour of the area above the transition very accurately. It also yields the correct asymptotic behaviour for the area in the limit of high pressures. In this approach, the constraint of fixed link length in treated exactly while the closure constraint is satisfied in a mean field sense. However, such a mean-field theory fails to describe the collapsed phase, also yielding incorrect results for the case of nonzero .
In Sec. V.2, we generalise an earlier mean field theory for the freely jointed chain to include semi-flexibility, imposing the constraint of fixed bond length via a Lagrange multiplier Gaspari et al. (1993); Haleva and Diamant (2006). The closure condition is imposed exactly. We thus derive expressions for the average area of the ring for all pressures and bending rigidity.
v.1 Density matrix mean-field for flexible polymers
In variational theory, a trial density matrix is chosen to approximate the actual density matrixChaikin and Lubensky (1995). The variational parameters are determined by minimising the variational free energy with respect to the parameters. The simplest mean-field theories assume a trial density matrix that is a product of independent single particle matrices, i.e,
where is the single particle density matrix of particle . The variational mean-field free energy is
The variational form for the density matrix should satisfy the constraint .
We choose the single particle density matrix based on the high pressure limit. In this limit, the ground state of our Hamiltonian is a regular N-gon, where the angle of the tangent vector is . The single particle density matrix has a delta function peak at this value. At intermediate pressures, we therefore take the form of the density matrix to be a gaussian of width (the variational parameter) centered about :
where the normalisation ensures that and is the error function defined as
Using this form of the density matrix, we obtain
When , the pressure and bending terms in Eq. (37) can be combined, and the problem is equivalent to one of a flexible polymer () with an effective pressure .
The variational parameter is chosen to be the that minimises in Eq. (37). This is done numerically. The average area, equal to , is then given by
We now derive the asymptotic behaviour of area in the limit of high pressures. We work in the limit when is large. For large pressures, we expect that tends to zero. In this limit
and the variational free energy is then given by
where . Solving , it is straightforward to obtain
The area then reduces to
For flexible polymers ( ), this mean-field theory reproduces the behaviour very accurately. It also obtains the correct asymptotic behaviour. In Fig. 6, we compare the Monte Carlo data for with the results of the above mean field theory and contrast it with the meanfield theory of Ref. Haleva and Diamant (2006).
The density matrix mean-field however, fails to correctly obtain the behaviour for non-zero values of the bending rigidity. It predicts a first order transition for , in disagreement with results from scaling theory. We compare the results of this mean field with the Monte Carlo data in Fig. 7 for a system with . This mean-field approach then predicts a transition at . The discrepancy between the two curves increases for larger values of .
We now describe an alternative mean-field approach to this problem which extends the harmonic spring-based mean field theory of Ref. Haleva and Diamant (2006) to non-zero values of .
v.2 Harmonic spring mean-field for semiflexible polymers
We follow the approach of Refs. Gaspari et al. (1993); Haleva and Diamant (2006) wherein the rigid links between particles are replaced by extensible springs. The spring constant of the springs is identified with a Lagrange multiplier, chosen so that the mean length of a spring equals unity.
Consider a partition function for N particles given by,
Note that while pressure couples to , the bending rigidity couples to the unit vectors . We make the approximation of replacing by . This makes the problem analytically tractable.
Expanding the tangent vectors in Fourier space as,
where , . The partition function then reduces to
By completing the squares, this integral can be written as a gaussian integral and hence can be calculated exactly. This gives
The parameter is determined by equating the mean square link length to one, i.e
When , the first factor in Eq. (49) becomes independent of , and then the resultant expression can be evaluated exactly. Hence, an analytic expression for can be obtained in this case Haleva and Diamant (2006). For , this is no longer possible, and for finite system sizes the resultant equation must be solved numerically. When , it is still possible to extract the behaviour of the system analytically.
We now determine the phase boundary from Eq. (49). We will consider the limit . First, note that for all . For positive , this gives the condition that . Second, consider the term in the denominator for . It is . If we assume that is continuous in , we have the second constraint that .
Setting and converting the first sum in Eq. (49) to an integral, the equation for reduces to
The critical pressure is obtained when the ratio becomes equal to , i.e.
For large values of , this goes as , which differs by a factor of from the answer obtained by scaling arguments [see Eq. (26)].
We shall now estimate in the different scaling regimes. We assume that is a non-decreasing function of (as in ). Then, since we have the constraint of , the ratio must continue to remain at for . Thus, above the critical point, we obtain
However, a simple substitution of Eq. (53) in Eq. (49) for does not satisfy Eq. (49). We therefore need to calculate the correction term arising from large but finite . We start by considering Eq. (49). The first term can be summed exactly, giving
We calculate the finite size corrections to as follows. Let
When , the main contribution to the left hand side of Eq. (54) comes from the term. The contribution from other is convergent as . Expanding the right hand side as a series in , we obtain
The independent term in the right hand side of Eq. (56) is nonzero for and is equal to zero for . Thus, when , we keep only the first term in the right side, while at , we need to keep the second term too. Solving for , we obtain,
We are now in a position to calculate the mean area from . This gives,
The numerical values obtained for are then substituted in this equation to get the corresponding value of the area. We can, however, analytically determine the scaling behaviour of the area in the limit of large system sizes from the values of calculated above.
For , we have,
This mean field theory reproduces the qualitative behaviour of the simulation data correctly. It predicts a continuous transition for all , unlike the density matrix field theory. However, there is a quantitative disagreement with the data. This can be seen by comparing the results of this mean-field theory with the simulation data in both the flexible (Fig. 6) and semi-flexible (Fig. 7) polymer cases.
Vi Scaling and Critical Exponents
The order parameter that describes the collapsed to inflated phase transition is the ratio of the area to the maximum area. When , the ratio is zero below the transition and non-zero above it. The behaviour near the transition line can be described by the scaling form
To obtain the one independent exponent, we resort to the scaling theory (see Sec. IV). At , . At the critical point, the area scales as . Combing with Eq. (63), we obtain and . These exponents are independent of .
We now look at the fluctuations. Consider the compressibility defined as
Thus, diverges as below the transition point. The behaviour near the transition point is described by the scaling form
where is a scaling function and . When , then . When , then . Comparison with Eq. (65) gives .
In Fig. 9, we plot the compressibility scaled as in Eq. (66) for two different values of . A good collapse is obtained again showing that the Flory type scaling theory gives the correct exponents. Similar, but noisier data can be obtained for the discrete model. We thus conclude that the introduction of semiflexibility does not affect any of the exponents describing the transition.
Vii The lattice problem
In this section, we present some additional enumeration results for the lattice problem. Consider the scaling theory presented in Sec. IV. The inextensibility of the polymer was captured by the term for a polymer of extent . This was obtained from a calculation based on the extension of a polymer under a force. Here we present numerical evidence supporting this.
where the scaling function is given by
We consider the corrections to the scaling form in Eq. (67). Let
Scaling theory predicts that should be a function of one variable . This is verified in Fig. 10 where is plotted against for a range of system sizes.
We also study the behaviour of area when is very large. When , the behaviour is seen to differ from the discrete version of the problem. It can be shown to beMitra et al. ()
In this paper, we have proposed and studied lattice and discrete models for self-intersecting pressurised semi-flexible polymers. Our work generalises results of Ref. Haleva and Diamant (2006) to include a bending rigidity. A simple variational mean-field approach provides very accurate fits to the Monte Carlo data for this problem in the absence of semi-flexibility. The mean-field approach for Rudnick and Gaspari (1991); Gaspari et al. (1993); Haleva and Diamant (2006) was generalised to the semiflexible case. The phase boundary between collapsed and inflated phases as well as expressions for the area as a function of and in the different phases were obtained analytically.
We have shown that the essence of the physics is captured through simple Flory approximations. The scaling predictions of the Flory theory were verified numerically for both the lattice and discrete cases.
We have also investigated the behaviour of the system in the extreme limits of a fully pressurised polymer ring and a collapsed configuration. For the fully pressurised ring, we deduce the leading order asymptotic behaviour of the area in both the discrete and lattice cases. The collapsed phase was studied by mapping this problem onto a quantum mechanical problem of an electron confined to two dimensions and placed in a transverse magnetic field. The analytic results thus obtained fit the data accurately.
The usefulness of these results for more realistic systems lies in the fact that both the restriction to the signed area as well as allowing for self-intersections at no energy cost are irrelevant in the large limit. The results obtained at large should therefore apply both qualitatively and quantitatively to the more realistic case of a pressurised self-avoiding polymer, where the pressure term couples to the true physical area and not to the signed area. This is the LSF model Leibler et al. (1987). The approach presented here is thus also useful in understanding the behaviour of a larger class of models, some of which are more physical in character, but which lack the analytic tractability of the model proposed and studied here.
Acknowledgements.This work was partially supported by grant 3504-2 of the Indo-French Centre for the Promotion of Advanced Research and the DST, India (GIM).
Appendix A Analytic answer in the small pressure regime
The problem of self-intersecting polymers in two dimensions with no bending rigidity () is analogous to the quantum mechanical problem of an electron moving in a magnetic field applied transverse to the plane of motion. Using this analogy, we obtain analytic expressions for the partition function and , the number of closed walks of area .
When an electron goes around a magnetic field, it picks up a phase factor proportional to the flux enclosed by the path. This flux is proportional to the product of the strength of the magnetic field times the algebraic area enclosed by the loop. The propagator then is the sum over all such loops. This suggests that the partition function for the polymer problem can be obtained from the quantum mechanical propagator for the electron problem provided the constants are appropriately mapped.
For an electron of charge and mass in a constant external magnetic field , in the z direction, in the case when the electron returns to the origin, the kernel can be written as Feynman and Hibbs (1965),
where . It picks up a flux given by
The motion of a quantum mechanical particle is governed by the two-dimensional Schroedinger equation,
The classical diffusion equation for a particle in two-dimensions is
where is the probability of finding the particle at . Thus to map the results of the quantum problem onto the polymer problem, we need to map and also identify