We study a solution of interacting semiflexible polymers with curvature energy in poor-solvent conditions on the d-dimensional cubic lattice using mean-field theory and Monte Carlo computer simulations. Building upon past studies on a single chain, we construct a field-theory representation of the system and solve it within a mean-field approximation supported by Monte Carlo simulations in d = 3. A gas-liquid transition is found in the temperature-density plane that is then interpreted in terms of real systems. Interestingly, we find this transition to be independent of the bending rigidity. Past classical Flory–Huggins and Flory mean-field results are shown to be particular cases of this more general framework. Perspectives in terms of guiding experimental results towards optimal conditions are also proposed.

Determining the phase behaviour of a solution of flexible and semiflexible polymers in poor-solvent conditions is a particularly challenging problem for several reasons. Unlike the case of colloidal liquids, where unambiguous gas-liquid and liquid-solid transitions are theoretically well characterized1 and experimentally observed,2 in the case of polymer solutions the presence of chain connectivity3–6 makes a full understanding of its phase behavior much more challenging, in particular at high concentrations.

One of the emerging conceptual problems hinges on the difficulty to discriminate between purely kinetic effects and those associated with the underlying thermodynamics.7,8 For instance, polymers do not completely crystallize when cooled down but become structured into a hierarchy of ordered structures.9 Also, it has been argued that on cooling a polymer melt undergoes a spinodal decomposition thus making the crystallization metastable and leaving the system out of equilibrium.7 Another difficulty stems from the large number of thermodynamic and structural parameters that need to be taken into consideration: in fact, in addition to usual thermodynamic quantities such as temperature and pressure that control the system, many other microscopic parameters such as interchain (in addition to the intrachain) interactions, the number of monomers in a chain, the stiffness of the fiber and the total polymer volume fraction have to be taken into account3–6 and become axes of a large parameter space.

Particularly important appears to be the case of semiflexible polymers, as a paradigmatic example of the protein folding problem10–13 or biopolymers in a crowded environment,14 and this is currently stimulating many studies along these lines. Our current understanding of these system comes, in particular, from computer simulations which, however, have been limited so far to particular situations. A first group of studies derived the gas-liquid phase diagram for flexible15 and semiflexible16 bead-spring chains up to only 100 monomers per chain using Monte Carlo (MC) simulations. In particular, they found a phase-equilibrium diagram very similar to that of simple liquids with some minor effects ascribed to the bending rigidity. Similar results were obtained by more recent and extensive simulations.17 Other simulations aimed at understanding entanglement properties between distinct chains18–20 or the onset of nucleation process.21 However, a comprehensive picture of the phase behavior of semiflexible polymers in poor-solvent solutions is currently lacking.

Surprisingly, even in the case of a single semiflexible polymer, a general theoretical understanding of the phase behaviour is still lacking notwithstanding several studies with different techniques have recently appeared22–27 that remained, however, focused on rather specific questions. For instance, by using mean-field arguments supported by bead-spring MC simulations, it has been argued that the ground state of a single semiflexible chain can be either a rod-like or a toroidal structure depending on the bending rigidity and the contour length of the polymer28,29 and this has been confirmed recently30 by computer simulations which also accounts for the temperature dependence. Interestingly, classical studies of a lattice model31,32 observe only rod-like phases, in the form of Hamiltonian paths, likely due to the geometrical constraints imposed by the lattice. Remarkably, the observed phase diagram is in excellent agreement with mean-field predictions using a field theoretical approach.13 

While a full understanding of the differences observed in lattice and off-lattice simulations is an interesting issue on its own right (and that will be discussed elsewhere), the present study will focus on providing the multi-chains counterpart of the aforementioned single chain studies.13,31,32 Specifically, in the wake of the long-standing tradition of lattice models for modeling polymer structure,10–13,33–45 here we develop a field-theoretical description of semiflexible self-avoiding chains with attractive interactions on the d-dimensional hypercubic lattice, and solve it within a mean-field approximation. Lattice grand canonical MC simulations will then be presented demonstrating the remarkable accuracy of the mean-field predictions. The multi-chain field theory approach that is proposed here build on past work by des Cloizeaux35 that extends the classical relation pointed out long ago by de Gennes4,34 of self-avoiding walks as the n → 0 limit of a spin O(n)-model, n being the number of components of each spin on a lattice point.

Notwithstanding its limitation, the present lattice approach has the great merit of making transparent the underlying physics and provide a guidance of the regions of this large parameter space that warrant a more in-depth analysis with dedicated numerical or experimental techniques. As a by-product of the theoretical analysis, we will re-obtain some classical results within a wider framework, as we will see.

The paper is organized as the following. In Sec. II, we provide a concise summary of the current understanding of the phase behavior of interacting semiflexible polymers. In particular, due to its specialized nature and for having inspired this work, we are going to highlight some salient aspects related to the past work of Orland and colleagues.10,13,42 Moreover, we will also discuss some conclusions from recent computational work by other groups which has the advantage of providing a broader view beyond the mean-field approach. The novel part of this work starts in Sec. III A, where we will introduce the exact grand canonical partition function Z of the lattice model for a multi-chain system (i.e., a polymer solution) on the d-dimensional cubic lattice which takes into account the local bending stiffness of the polymer fiber, excluded-volume and short-range attractive interactions between close-by monomers. Then, by exploiting the analogy between self-avoiding polymers and the spin O(n → 0)-model, we construct (Sec. III B) the exact field-theoretic representation of Z. Since an exact computation of Z is unfeasible, we describe a uniform saddle-point approximation (Sec. IV) and obtain the corresponding mean-field solution of the problem, the reliability of which is demonstrated by comparison to MC computer simulations (Secs. V and VI). Finally we show that our results (Sec. VII) recapitulate, as particular cases, several models that have been discussed in the past and we demonstrate the equivalence between our approach and the classical Flory–Huggins theory46,47 for mixtures. Discussion and conclusions, with an outline on open problems and possible future perspective, are presented in Sec. VIII.

Before introducing (Sec. III) our field-theoretic formalism for semiflexible polymer solutions and in order to set the stage, it proves instructive to recapitulate the single chain formalism.10,13,42

Denoting by L the linear size of the hypercubic lattice in d dimensions and by a the lattice spacing, we consider Hamiltonian paths, polymer chains whose number of nodes N is equal to V/ad = (L/a)d where V is the volume of the lattice: that is, we consider all lattice points occupied and no vacancies. For simplicity, we further restrict our considerations to closed paths, or Hamiltonian rings (HR), knowing that the statistics is the same in the thermodynamic limit.42 

Consider the following O(n)-vector model where the n-component vector,
S(x)S1(x),S2(x),,Sn(x),
(1)
is associated to each lattice point x. By defining the scalar product S(x)S(x)i=1nSi(x)Si(x) between any two vectors associated to lattice points x and x, we assume the following constraint on the norm-square of S(x):
S(x)2S(x)S(x)=i=1nSi(x)2=n.
(2)
The reduced Hamiltonian for the n-vector model in zero external field reads
βH=J2x,xΔ(x,x)S(x)S(x),
(3)
where β = 1/(kBT), kB being the Boltzmann constant and T being the temperature, J is the coupling constant between the spins, and where
Δ(x,x)=1,if|xx|=a,0,otherwise.
(4)
From now on we assume periodic boundary conditions and Eq. (3) can be also written as
βH=Jxμ=1dS(x)S(x+êμ),
(5)
where êμ is the unit lattice vector pointing towards the “+μ”-direction. Then, we introduce the integration measure,
dΩn(x)=dS(x)δi=1nSi(x)2n,
(6)
with the Dirac δ-function enforcing the constraint Eq. (2) and the related geometrical average (or, trace operation)
ΩxdΩn(x)()xdΩn(x).
(7)
A very peculiar feature of the operation defined in Eq. (7) is that, in the formal limit n → 0, the following equality for the moment-generating function holds:4 
xeS(x)φ(x)Ω=x1+12φ(x)2,
(8)
i.e., the moment-generating function has a simple quadratic form. Further consequences of this fact will be illustrated in  Appendix A.
Let us consider now the quantity eβHΩ. One can show that, in the formal limit n → 0, the following expansion holds:
eβHΩ=1+n=1NJZ+On2,
(9)
where Z is the total number of self-avoiding closed paths of total length [extended details on the derivation of Eq. (9), which is non-trivial, are provided in  Appendix A]. At the same time, by using the standard Hubbard–Stratonovich transformation,48,49 we obtain
eβHΩ=xdφ(x)eAeJxS(x)φ(x)Ωxdφ(x)eA,
(10)
where
A=12x,xΔ1x,xφ(x)φ(x),
(11)
with Δ−1 the inverse of the matrix Δ [Eq. (4)]. At this point it has to be noted that, strictly speaking, the matrix Δ is not positive definite, therefore the Hubbard–Stratonovich transformation itself is, in principle, ill-defined. However, this technical difficulty can be overcome through the more rigorous approach42,44 involving Fresnel integrals which leaves the final results unaffected. In the end then and in the limit n → 0, the numerator of Eq. (10) can be easily computed by resorting to Eq. (8), granting the result
eβHΩ=xdφ(x)eAx1+J2φ(x)2xdφ(x)eA.
(12)
We now notice that ZN in Eq. (9) – which coincides with the total number of HR on our lattice – can be formally50 obtained as
ZN=limn0limJ1n1JNeβHΩ.
(13)
Finally, by combining Eq. (12) with (13) we get the following expression for ZN,
ZN=limn01nxdφ(x)eAx12φ(x)2xdφ(x)eA,
(14)
which was introduced first in Ref. 42. It is important to stress that, in order to compute ZN, we have used the fact that the trace operation Eq. (7) has the very peculiar properties described in  Appendix A. Finally, for the purpose of computing ZN, the geometrical origin of this trace (i.e., constraining the spin vectors on the surface of a sphere of radius n) becomes completely irrelevant.

We demonstrate now that there exists an alternative method51 of finding ZN (14) that represents both a shortcut with respect to the approach presented so far and has the great advantage of being exportable to the more general situation considered here (Sec. III) in a relatively straightforward manner.

The method consists in defining a priori a trace operation – that we denote by the symbol ⟨·⟩0 – characterized by the desired mathematical properties:
10=0,
(15)
Si0=0,
(16)
SiSj0=δij,
(17)
Si1Si2Sik0=0,ifk3,
(18)
between spin components on the same lattice site, while S-vectors on different sites are independent from each other. Notice that the only difference from the trace (7) is that now the trace of 1 is equal to 0 (see  Appendix A).
Based on the definitions (15)(18) and by taking J = 1 in the Hamiltonian of the n-vector model [Eq. (3)], the partition function (14) of the HR is equivalent to:
ZN=limn01neβH0.
(19)
In fact, by taking the Hubbard–Stratonovich transformation of the term inside brackets [Eq. (10) with J = 1], we have
ZN=limn01nxdφ(x)eAxeS(x)φ(x)0xdφ(x)eA=limn01nxdφ(x)eAx1+(S(x)φ(x))+12(S(x)φ(x))20xdφ(x)eA=limn01nxdφ(x)eAx12φ(x)2xdφ(x)eA,
(20)
and the last line of Eq. (20) is identical to Eq. (14). Notice that the factorization of the product in the second line of Eq. (20) and the expansion up to second-order follow straightforwardly from definitions (15)(18). This concludes the proof.

In Sec. III, we will present a suitable generalization of definitions (15)(18) to treat solutions of semi-flexible polymers with bending stiffness and monomer-monomer attractive interactions for poor-solvent conditions.

We generalize here the formalism introduced in Sec. II and we consider a system (i.e., a solution) of semiflexible linear polymer chains with attractive interactions between non-bonded monomer pairs modelling poor-solvent conditions.6 

Chains are arranged on the same hypercubic lattice in d dimensions introduced in Sec. II. Again, lattice spacing, linear side length, volume and total number of sites are denoted, respectively, by: a, L, V = Ld and N = V/ad. Chains are self- and mutually-avoiding, i.e., any two monomers – be they from the same or different chains – cannot occupy the same lattice site. Chain stiffness is modeled by introducing a bending energy penalty ϵa > 0 for two consecutive bonds along the same chain forming a turn (or, an angle), while attractive interactions between non-bonded monomers are accounted for by an energy reward −ϵi < 0 (ϵa, ϵi > 0) for any two monomers which are separated by a unit lattice distance and either are non-consecutive if they belong to the same chain or they are on distinct chains.

For computational convenience, we work in the grand canonical ensemble where neither the number of chains nor the number of bonds are fixed and we introduce the grand canonical partition function,
Z=Z(κ,η,ϵa,ϵi)={C}w(κ,η,ϵa,ϵi;C),
(21)
where the sum is intended over the set of all possible configurations {C} and where the thermal (Boltzmann) weight for each conformation C,
w(κ,η,ϵa,ϵi;C)=κNb(C)η2Nc(C)eβϵaNa(C)+βϵiNi(C),
(22)
depends on: (a) Nb(C), the total number of bonded monomer pairs with corresponding bond fugacity κ; (b) Nc(C), the total number of chains with corresponding chain fugacity52  η2; (c) Na(C) and Ni(C), respectively the total number of corners and the total number of non-bonded monomer pairs separated by one lattice distance. Again, β = 1/(kBT) is the Boltzmann factor at temperature T and kB is the Boltzmann constant, and we assume periodic boundary conditions through the hypercubic lattice. Less obviously, we anticipate here and justify briefly in Sec. III B that closed chains are ruled out in our field theory. An example illustrating a particular configuration C on the square lattice (d = 2) is shown in Fig. 1. Notice that, per our definition, the smallest length of a single chain corresponds to 1 lattice bond.
FIG. 1.

Schematic illustration of a particular configuration C on the square lattice (d = 2). By assuming periodic boundary conditions (PBC, see text for details), we have: Nc(C)=4 chains, Nb(C)=27 bonds, Na(C)=14 turns or angles (marked as blue corners), and Ni(C)=15 pairs of interacting monomers (connected by dashed red lines, with two pairs interacting through PBC.

FIG. 1.

Schematic illustration of a particular configuration C on the square lattice (d = 2). By assuming periodic boundary conditions (PBC, see text for details), we have: Nc(C)=4 chains, Nb(C)=27 bonds, Na(C)=14 turns or angles (marked as blue corners), and Ni(C)=15 pairs of interacting monomers (connected by dashed red lines, with two pairs interacting through PBC.

Close modal

The central quantity of our work, the grand canonical partition function Z (21), admits a field-theoretic representation.

To show it, the first point consists in devising a method “to count” the total number of bonds (Nb), chains (Nc), angles (Na), and non-bonded monomer-monomer pairs (Ni) characterizing each given chain configuration C. To this purpose, we start by defining the scalar function for the configuration C,
ωC(x)=1,if lattice positionxis occupied by a monomer0,otherwise
(23)
By using the matrix Δ(x,x) [Eq. (4)], we have
Nb(C)+Ni(C)=12x,xΔ(x,x)ωC(x)ωC(x),
(24)
and the Hubbard–Stratonovich transformation48,49 of the exponential of the rhs of Eq. (24) is equivalent to the expression containing the scalar field ψ=ψ(x)
Dψexp12x,xΔ1(x,x)ψ(x)ψ(x)+xωC(x)ψ(x)
(25)
with
Dψ(2π)N/2(detΔ)1/2xdψ(x).
(26)
Then, at each lattice position x we introduce d distinct n-component vectors, Sμ(x)(Sμ1(x),Sμ2(x),,Sμn(x)) with μ = 1, 2, …, d, obeying the generalized trace rules:
10=1,
(27)
Sμi0=0,
(28)
SμiSνj0=δij[δμν+(1δμν)eβϵa],
(29)
Sμ1i1Sμ2i2Sμkik0=0,ifk3.
(30)
Again, S-vectors on different sites are independent from each other under the trace operations just defined.
By using Eqs. (24)(26) and the discussion in  Appendix A,
Z=Dψe12x,xΔ1(x,x)ψ(x)ψ(x)×limn0x1+H(x)μ=1dSμ1(x)×μ=1d1+h(x)h(x+êμ)Sμ(x)Sμ(x+êμ)0,
(31)
where êμ is the unit lattice vector pointing towards the “+μ”-direction introduced in Sec. II and
H(x)=η1+(d1)eβϵaeβϵi2ψ(x),
h(x)=κeβϵi/2eβϵi2ψ(x).
Importantly, it must be noticed that Eq. (31) takes into account the fact that there must be no branching points [owing to the trace definitions (27)(30), any branching point gives a contribution equal to 0] and no closed loops (from  Appendix A, the statistical weight of any configuration C with k closed-loops is proportional to nk and so its contribution disappears in the n → 0 limit). In fact, the term of Eq. (31) appearing under the limit of n → 0 can be also written as
limn0x1+H(x)μ=1dSμ1(x)×expμ=1dh(x)h(x+êμ)Sμ(x)Sμ(x+êμ)0=limn0x1+H(x)μ=1dSμ1(x)×exp12x,xμ=1dΔμ(x,x)h(x)h(x)Sμ(x)Sμ(x)0,
(32)
since, again because of the trace definitions (27)(30), in the expansion of the exponential all terms higher than the first are of order n and they disappear in the n → 0 limit and where, by analogy with the matrix Δ [Eq. (4)], we have defined
Δμ(x,x)=1,if|(xx)êμ|=a,0,otherwise.
(33)
For expression (31) to be helpful, we have to remove the dependence on the S-vectors in favor of real fields. To this purpose, we perform a Hubbard–Stratonovich transformation48,49 of the exponential term of the last line of Eq. (32) containing the S-vectors,
e12x,xμ=1dΔμ(x,x)h(x)h(x)Sμ(x)Sμ(x)=Dφe12x,xμ=1dΔμ1(x,x)φμ(x)φμ(x)+xh(x)μ=1dSμ(x)φμ(x),
(34)
where we have introduced the d real vector fields φμ(x) (μ = 1, …, d) with φμ(x)(φμ1(x),φμ2(x),,φμn(x)) and the corresponding measure [see Eq. (26), for analogy]
Dφ(2π)ndN/2μ=1ddetΔμn/2xμ=1ddφμ(x).
(35)
Finally by (i) inserting Eq. (34) into Eq. (31) via Eq. (32); (ii) Taylor-expanding the term containing the φμ-fields; (iii) applying the trace definitions (27)(30) and (iv) noticing that the first two terms in Eq. (35) give = 1 in the limit n → 0, one can show that, up to an unimportant multiplicative constant,
Z=limn0xdψ(x)xμ=1ddφμ(x)×expA[{ψ}]μ=1dAμ[{φμ}]+xln1+eβϵiψ(x)B[{φμ(x)}],
(36)
where we have defined the following functionals:
A[{ψ}]=12x,xΔ1(x,x)ψ(x)ψ(x),
(37)
Aμ[{φμ}]=12x,xΔμ1(x,x)φμ(x)φμ(x),
(38)
B[{φμ(x)}]=κeβϵi2(1eβϵa)μ=1d|φμ(x)|2+eβϵaμ=1dφμ(x)2+κηeβϵi/2μ=1dφμ1(x).
(39)
Importantly, notice the explicit presence of φμ1(x) in Eq. (39). This is a direct consequence of the fact that, in order to describe a system of multiple chains through the O(n → 0) formalism, it suffices to introduce an external magnetic field in the spin Hamiltonian.4,35 This field can pick any arbitrary direction: in our derivation, we have chosen the direction with μ = 1. As a validation of Eqs. (36)(39), we report that, in the “single-chain” limit η → 0, we get back the original result by Doniach et al.13 for a single semiflexible chain with non-bonded attractive interactions in the presence of lattice vacancies.

The exact grand canonical partition function Z [Eq. (36)] is the central result of this work. As for the single chain case a direct evaluation of Z is not feasible but the field theoretical formulation [Eq. (36)] is very suitable for its mean field (MF) treatment.10,13,42

We start by differentiating the exponential in Eq. (36) with respect to φμi(x) and ψ(x) and set the obtained expressions equal to 0 in order to get the stationary solution. We further take the solutions to be homogeneous assuming translational invariance and break the O(n) symmetry of the vector field so that
φμ(x)=(φ,0,,0),
(40)
ψ(x)=ψ,
(41)
for every x and every μ, thus obtaining
φ2=eβϵiψκeβϵiq(β)2φ+κηeβϵi/21+eβϵiψdκeβϵiq(β)4φ2+κηeβϵi/2φ,
(42)
ψ2d=βϵieβϵiψdκeβϵiq(β)4φ2+κηeβϵi/2φ1+eβϵiψdκeβϵiq(β)4φ2+κηeβϵi/2φ,
(43)
where13, q(β)=2+2(d1)eβϵa. In terms of the solutions53  φ = φ(κ, η, ϵi, ϵa) and ψ = ψ(κ, η, ϵi, ϵa) of the MF Eqs. (42) and (43), the grand potential per lattice site (up to unimportant additive constants) reads
βΩ(κ,η,ϵa,ϵi)=ψ24d+dφ24ln1+deβϵiψ×κeβϵiq(β)4φ2+κηeβϵi/2φ.
(44)
Notice that with the ansatz (40) and (41) every dependence upon n disappears, and thus the limit n → 0 is trivial.

On setting η = 0, Eqs. (42) and (43) reduce to the ones obtained in Ref. 13 for the single chain model. In the following, we will thus solve the saddle-point equations in the case η > 0 that will be then compared with Monte Carlo simulations (Sec. V) in Sec. VI.

In order to check the validity of the MF approximation as well as to assess its limits, we have performed Metropolis54 grand canonical Monte Carlo (GCMC) computer simulations of the lattice model (Sec. III) on the three-dimensional cubic lattice. Essentially, the goal of the GCMC simulations is to obtain a representative sample of polymer configurations in agreement with the grand canonical partition function (21).

The implementation of our algorithm is relatively straightforward, and it works as the following. As explained in Sec. III, the Boltzmann weight w [see Eq. (22)] of each polymer configuration C in the ensemble is a function of the total number of bonds [Nb(C)], distinct chains [Nc(C)], turns [Na(C)] and pairs of non-bonded nearest-neighbor monomers [Ni(C)]. Therefore, at each MC step one single bond is randomly inserted in or removed from the lattice, provoking a change of the configuration C0 to the configuration C1. In order to enforce the condition of detailed balance, we accept54 the new conformation with probability given by the expression:
acc(C0C1)=min1,dϕb,1w(κ,η,ϵa,ϵi;C1)w(κ,η,ϵa,ϵi;C0)(bond inserted),min1,ϕb,0dw(κ,η,ϵa,ϵi;C1)w(κ,η,ϵa,ϵi;C0)(bond removed),
(45)
where ϕb,0 (respectively, ϕb,1) is the bond density [see Eq. (46)] of the configuration C0 (resp., C1). Whenever the insertion of a new bond leads to a forbidden configuration (e.g., for the presence of branching points or closed loops), the move is automatically discarded.

In order to check for finite-size effects, we have performed preliminary calculations and compared corresponding results for lattice sizes L/a = 4, 8, 16. This analysis indicates that for L/a ≥ 8 the results do not change significantly with the lattice size and hence we will fix L/a = 8 in all calculations henceforth. This guarantees a good compromise between computational efficiency and accuracy.

For each chosen pair of (κ, η) values, we have led a simulation run consisting of 107 MC steps. Then, for each GCMC trajectory, a standard block analysis procedure55 has been carried out in order to estimate uncertainties on the considered physical observables. Every trajectory has been checked individually in order to make sure that all the curves obtained from corresponding block analyses and representing the MC-time evolution of the distinct quantities have completely equilibrated. This procedure has been applied to the number of bonds, the number of chains and the internal energy of the system.

From our mean-field estimate of the grand potential, Eq. (44), we can compute:

  • The bond density,
    ϕbNbN=βκΩ(κ,η,ϵi,ϵa)κ=d4φ2;
    (46)
  • The chain density,
    ϕcNcN=βη2Ω(κ,η,ϵi,ϵa)η=d2eβϵiψκηeβϵi/2φ1+deβϵiψκeβϵiq(β)4φ2+κηeβϵi/2φ;
    (47)
  • The total monomer density,
    ϕϕb+ϕc=1βϵiψ2d.
    (48)
    In the rest of this section, we specialize the saddle-point Eqs. (42) and (43) to various particular cases and compare the corresponding results to Monte Carlo simulations in d = 3.

The simplest case to be considered is the case of non-interacting (still non-overlapping) flexible chains with no bending penalty nor monomer-monomer attractive interactions. In spite of its simplicity, this case proves to be rather instructive. In this case Eqs. (42) and (43) read
φ2=κdφ+κη1+dκd2φ2+κηφ,
(49)
ψ2d=0.
(50)
The only relevant field is thus φ and, since it satisfies the simple cubic Eq. (49), in principle three (real) solutions can be possible. However two additional constraints identify the only acceptable solution. First, the argument of the logarithm in Eq. (44) [i.e., the denominator in Eq. (49)] must be strictly >0. Second, the chain density [Eq. (47)] must satisfy the inequality 0 ≤ ϕc ≤ 1/2. In  Appendix B we provide evidence that, for every κ ≥ 0 and η > 0, there exists one and only one such solution φ > 0 which is a continuous function of the parameters. This solution can then be inserted in Eqs. (46) and (47) to obtain the bond and chain density, ϕb and ϕc.

MF calculations for the bond and chain density, ϕb and ϕc, as a function of the bond fugacity κ and for two representative chain fugacities η = 0.2 (small) and η = 1.5 (large) are shown as solid lines in Fig. 2 [panels (a) and (b), respectively] and compared to corresponding GCMC simulations (symbols). The nearly quantitative agreement between the MF calculations and the GCMC simulations is remarkable, thus validating our MF approach.

FIG. 2.

ϵa = ϵi = 0. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are, respectively, for the MF solution and the GCMC computer simulations.

FIG. 2.

ϵa = ϵi = 0. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are, respectively, for the MF solution and the GCMC computer simulations.

Close modal

One striking feature of the bond density curves [see Fig. 2(a)] is that they intersect at a certain κ = κ*, such that ϕb(κ*) ≃ 0.5. Although odd at first sight, this behaviour can be simply rationalized as the following. When ϕb < 0.5 it is likely that the insertion of a new bond will also lead to the creation of a new chain. Thus, for κ < κ*, the bond density increases faster for larger values of η [red curve in Fig. 2(a)] than for smaller values of η (blue curve) because configurations with a larger number of chains are more favoured. Conversely, when ϕb > 0.5 once a new bond is inserted it will link two different chains, thus reducing their total number. Under these conditions, for κ > κ* the bond density increases faster for smaller values of η (blue curve) than for larger ones (red curve). A further support to this interpretation also stems by the fact that the chain density [Fig. 2(b)] has a maximum at κκ*.

It also proves instructive to derive simple analytical expressions for ϕb and ϕc in the limit κ ≫ 1. In this regime, Eq. (49) depends only on one parameter, namely the ratio η/κ, and has the following solutions
φ4d12dηκ,ifηκ12d+d8κη,ifηκ1
(51)
By plugging these results into the expressions for ϕb [Eq. (46)] and ϕc [Eq. (47)] we find
ϕb112dηκ,ifηκ112+d8κη,ifηκ1
(52)
and
ϕc12dηκ,ifηκ112d8κη,ifηκ1
(53)
For κ ≫ 1 the monomer density ϕ = ϕb + ϕc is always ≃1. However, there are two different scenarios: if η/κ1 then ϕb/ϕc ≫ 1, i.e., the number of chains is very low, but on average they are very long. On the other hand, if η/κ1 then ϕb/ϕc ≃ 1, that is the number of chains is large but they are all essentially formed by one single bond, in agreement with previous interpretation.
With only the contribution of the bending penalty but no monomer-monomer attractive interactions, Eqs. (42) and (43) read
φ2=κq(β)2φ+κη1+dκq(β)4φ2+κηφ,
(54)
ψ2d=0.
(55)
It is easy to see that this is the same situation of Sec. VI A with renormalized fugacities κκq(β)/2d and ηη2dq(β). Since 2 ≤ q(β) ≤ 2d, introducing a non-zero bending stiffness leads ultimately to a lower effective bond fugacity and a larger effective chain fugacity. Hence, as seen in Sec. VI A, again we have only one acceptable solution which is a continuous function of the parameters.

By fixing the bending stiffness to the convenient value ϵa = kBT, MF calculations for the bond and chain density, ϕb and ϕc, as a function of the bond fugacity κ and for chain fugacities η = 0.2 and η = 1.5 are shown as solid lines in Fig. 3 [panels (a) and (b), respectively] and compared to corresponding GCMC simulations (symbols). As in previous Sec. VI A, the agreement between MF calculation and GCMC simulations is remarkable.

FIG. 3.

ϵa/kBT = 1, ϵi = 0. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are as in Fig. 2.

FIG. 3.

ϵa/kBT = 1, ϵi = 0. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are as in Fig. 2.

Close modal

In principle, here one would have expected an isotropic-to-nematic transition which would be observable as a singularity in the average “angle” density, ϕa. However, due to the fact that the only contribution of the bending stiffness is to renormalize the fugacities, our MF treatment does not display the appearance of such transition. This is made evident in Fig. 4 where we show the average angle density as obtained in the MF approximation [solid lines, see formula (64)] in favorable comparison to the results of GCMC simulations (symbols).

FIG. 4.

ϵa/kBT = 1, ϵi = 0. Angle density ϕa as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are as in Fig. 2.

FIG. 4.

ϵa/kBT = 1, ϵi = 0. Angle density ϕa as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are as in Fig. 2.

Close modal
This is the complementary case of previous one, where monomer-monomer attraction is present but there is no bending penalty. In this case, Eqs. (42) and (43) read
φ2=eβϵiψ(κeβϵidφ+κηeβϵi/2)1+eβϵiψdκeβϵid2φ2+κηeβϵi/2φ,
(56)
ψ2d=βϵieβϵiψdκeβϵid2φ2+κηeβϵi/2φ1+eβϵiψdκeβϵid2φ2+κηeβϵi/2φ.
(57)
which do not admit a simple closed solution for ψ as in the previous cases. Hence these equation have to be solved numerically with the two constraints discussed in Sec. VI A plus the condition 0 ≤ ϕ ≤ 1, that implies 0ψ2dβϵi, for the total monomer density ϕ [Eq. (48)]. Interestingly, in this case one finds that there are multiple acceptable solutions for this system, and hence the most stable solution corresponds to that minimizing the grand potential Ω [Eq. (44)]. In particular, this leads to the appearance of discontinuities in ϕb and ϕc (see  Appendix C for details).

In fact, by fixing the monomer-monomer interaction to the convenient value ϵi = kBT, MF calculations for the bond and chain density, ϕb and ϕc, as a function of the bond fugacity κ and for chain fugacities η = 0.15 and η = 0.25 are shown as solid lines in Fig. 5 [panels (a) and (b), respectively] and compared to corresponding GCMC simulations (symbols). In both panels one can easily distinguish two different phases, with both densities acting as order parameters. In the first (gas-like) phase, the total monomer density ϕ = ϕb + ϕc is close to 0 and ϕb/ϕc ≃ 1. This is valid up to some critical value κcr above which ϕ is close to 1 and ϕb/ϕc > 1 (liquid-like phase).

FIG. 5.

ϵa = 0, ϵi/kBT = 1. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.15 (blue) and η = 0.25 (red). The discontinuity predicted by MF calculations (lines) is confirmed by GCMC simulations (symbols), yet the critical value of κ is slightly different between theory and simulations.

FIG. 5.

ϵa = 0, ϵi/kBT = 1. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.15 (blue) and η = 0.25 (red). The discontinuity predicted by MF calculations (lines) is confirmed by GCMC simulations (symbols), yet the critical value of κ is slightly different between theory and simulations.

Close modal

By varying systematically the parameters βϵi and η, we have extracted each corresponding critical value κcr through the numerical solution of the coupled Eqs. (56) and (57). Three illustrative “coexistence” lines corresponding to the values βϵi = 0.75, 1.00 and 1.25 are shown in panel (a) of Fig. 7 (solid lines). Similarly [panel (b), filled symbols], we have determined the values of the total monomer density ϕ [Eq. (48)] at the coexistence by varying βϵi systematically and for the two representative fugacity values η = 0.1 and η = 0.3. We will discuss these results in full fledged way in Sec. VIII.

When both the attractive monomer-monomer interaction and bending stiffness are non-zero, we need to solve the complete system of Eqs. (42) and (43). Again, the strategy is the same as in Sec. VI B: renormalizing the fugacities [κκq(β)/2d and ηη2dq(β)] to absorb the terms accounting for the bending stiffness and obtain a system of equations which is equivalent to that presented in Sec. VI C. The general behavior of the solutions will thus be exactly the same (and, so, the discontinuities in ϕb and ϕc) as for a system with no bending stiffness, the only difference being in the changing of the coexistence line between the two phases in the (κ, η)-plane, depending on the value of ϵa.

By fixing again the bending stiffness and the monomer-monomer interaction to the values ϵa = ϵi = kBT, MF calculations for the bond and chain density, ϕb and ϕc, as a function of the bond fugacity κ and for chain fugacities η = 0.15 and η = 0.25 are shown as solid lines in Fig. 6 [panels (a) and (b), respectively] and compared to corresponding GCMC simulations (symbols). As expected, GCMC simulations confirm MF calculations and we can distinguish once again the gas (ϕ = ϕb + ϕc ≈ 0) and liquid (ϕ ≈ 1) phases. Finally, analogously to Sec. VI C, we produce examples of coexistence and gas-liquid transition lines, see panel (a) [dashed lines, different colors are for different ϵi (see caption)] and panel (b) [empty symbols, different colors are for different η (see caption)] of Fig. 7.

FIG. 6.

ϵa/kBT = ϵi/kBT = 1. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.15 (blue) and η = 0.25 (red). As in the case with no bending stiffness (Fig. 5), the discontinuity predicted by MF calculations (lines) is apparent and its presence is confirmed by GCMC simulations (symbols). Again, the critical value for κ is quantitatively different between theory and simulations.

FIG. 6.

ϵa/kBT = ϵi/kBT = 1. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.15 (blue) and η = 0.25 (red). As in the case with no bending stiffness (Fig. 5), the discontinuity predicted by MF calculations (lines) is apparent and its presence is confirmed by GCMC simulations (symbols). Again, the critical value for κ is quantitatively different between theory and simulations.

Close modal
FIG. 7.

(a) Coexistence lines between the gas phase and the liquid phase in the (κ, η)-plane. Solid lines correspond to ϵa = 0, whereas dashed lines are for ϵa/ϵi = 1. Lines colors blue, green and red are for ϵi/kBT = 0.75, 1.00 and 1.25, respectively. Below the coexistence line the system is in the gas phase (ϕ ≃ 0), whereas above the coexistence line it is found in the liquid phase (ϕ ≃ 1). (b) Gas-liquid transition in the [ϕ, 1/(βϵi)]-plane. Above the critical point, the system is in a single homogeneous gas phase. Below the critical point, the system phase-separate in a gas phase coexisting with a liquid phase, and the figure display the coexistence (binodal) line. Symbols colors blue and red are for η = 0.1 and η = 0.3, respectively. Filled symbols correspond to ϵa = 0, whereas empty symbols are for ϵa/ϵi = 1. For each data set the relative grey symbol marks the value of the corresponding “critical” temperature (letters “C” and “C′”), with grey lines used for guiding the eye.

FIG. 7.

(a) Coexistence lines between the gas phase and the liquid phase in the (κ, η)-plane. Solid lines correspond to ϵa = 0, whereas dashed lines are for ϵa/ϵi = 1. Lines colors blue, green and red are for ϵi/kBT = 0.75, 1.00 and 1.25, respectively. Below the coexistence line the system is in the gas phase (ϕ ≃ 0), whereas above the coexistence line it is found in the liquid phase (ϕ ≃ 1). (b) Gas-liquid transition in the [ϕ, 1/(βϵi)]-plane. Above the critical point, the system is in a single homogeneous gas phase. Below the critical point, the system phase-separate in a gas phase coexisting with a liquid phase, and the figure display the coexistence (binodal) line. Symbols colors blue and red are for η = 0.1 and η = 0.3, respectively. Filled symbols correspond to ϵa = 0, whereas empty symbols are for ϵa/ϵi = 1. For each data set the relative grey symbol marks the value of the corresponding “critical” temperature (letters “C” and “C′”), with grey lines used for guiding the eye.

Close modal

It turns out that it is possible to obtain some information on the system even without directly solving Eqs. (42) and (43). As we are interested in the free energy rather than in the grand potential, we perform a Legendre transform in order to have a dependence upon ϕc and ϕb. To this aim, we need to be able to express κ and η in terms of these densities.

Making use of Eqs. (46) and (47), and without an explicit derivation of the solutions for φ and ψ, one can express the fugacities η and κ in terms of ϕb and ϕc as:
η=q(β)ϕcedβϵi(ϕb+ϕc)d(1ϕbϕc)(ϕbϕc),
(58)
κ=(ϕbϕc)eβϵiq(β)ϕb(1ϕbϕc)e2dβϵi(ϕb+ϕc).
(59)
This allows us to compute the reduced free energy per site
βf=βΩ+ϕblnκ+2ϕclnη
(60)
where it is clear that f depends on both, ϕb and ϕc. However, it is more convenient to express it in terms of ϕ and ̄=ϕ/ϕc (i.e., a measure of the average number of monomers per chain):
βf=dβϵiϕ2+(1ϕ)ln(1ϕ)+ϕ1̄lnϕ12̄lnq(β)e+11̄βϵi+12̄ln̄2̄11̄lnd̄(̄1)e,
(61)
where e = 2.718 28… is the Euler’s number. Equation (61) is a key result of our mean-field analysis. Notice that, since in our model the minimum length of a chain is defined to be = 1 bond (i.e., 2 monomers), it must be ̄2.
Now that we have a mean-field estimate of the free energy, we employ Eq. (61) and compute the internal energy U of the system:
UN=βfβ=NiN(ϵi)+NaNϵa.
(62)
where
NiN=dϕ2ϕ11̄,
(63)
NaN=q(β)2q(β)12̄ϕ.
(64)
Equations (63) and (64) bear interesting physical interpretations. Let us first discuss Eq. (63). Within our mean-field approach the number of interactions per lattice site NiN is not guaranteed to be always non-negative. From Eq. (63) in fact, we readily see this to be true only when ϕ1d11̄. On the other hand, in the limiting case ϕ = 1, it is easy to check that the density of interactions is exactly that predicted by Eq. (63), namely d(11/̄). Also, the interaction density is a decreasing function of the average chain length, as it should be since the interaction is only between non-consecutive nearest neighbour monomers. Then we focus on Eq. (64), noticing that the density of angles is linear in ϕ and the value of the bending rigidity can only modify the corresponding proportionality constant. From the dependence of Eq. (64) on the number of monomers per chain ̄, we see that: (i) if ̄=2 (i.e., all chains consist of 1 bond only) the total number of angles is = 0, as expected; (ii) increasing ̄ leads to a larger number of angles; (iii) the proportionality constant between the angle density and ϕ is less than 1 for every value of βϵa and ̄, as it should.
Further insights can be obtained by comparing our results with past work by Flory who also derived a statistical thermodynamics theory for semiflexible chain molecules on a lattice using combinatorial arguments.33 A central quantity of his theory is the mean bending degree of the chain,
gF=2(d1)eβϵa1+2(d1)eβϵa,
(65)
which results to be independent of the concentration. The same quantity can be computed [by using Eq. (64)] within our approach as well, i.e.,
g=NaN1ϕc(̄2)=(d1)eβϵa1+(d1)eβϵa,
(66)
is the average monomer fraction where the polymer chain displays a turn. As a direct consequence of Eq. (64), our estimate of g is also independent of the concentration ϕ, in agreement with Flory. This means that the degree of bending is only dependent on the temperature, irrespective of whether the chain is in a melt or in a dilute solution.56 

Note that a comparison between the two estimates gives g < gF. On the other hand we can also argue that gF is clearly overestimating the true value because it means that at each position along the chain, the possible directions available to make a turn are always 2(d − 1), thereby not accounting for the long-range correlations due to self-avoidance. By contrast, Eq. (66) implies that the number of possible directions to make a turn at each step along the chain is on average d − 1, thus effectively accounting for possible long-range correlations induced by self-avoidance.

As both estimates Eqs. (65) and (66) are independent of the concentration as remarked, it seems not plausible that it may act as an order parameter for a phase transition. Nevertheless, Flory postulates the existence of a phase-separated, ordered state (g = 0) also at T > 0, and assumes that the entropy of such state is 0. He then proceeds in computing the free energy difference between the latter and a completely mixed, disordered state (g ≠ 0), the free energy of which is computed within the same mean-field theory. He then verifies a posteriori that there exists a critical temperature below which the ordered state is thermodinamically more stable. The premise that the entropy of the ordered state is 0 is crucial for the derivation of Flory, and it is ultimately this assumption that leads to the appearance of a critical temperature, corresponding to a first-order phase transition. Later studies by Gujrati and co-workers36,38,40 were however able to derive an exact lower bound for the entropy of a system of self-avoiding chains on a lattice that was found to be strictly positive at any temperature T > 0, therefore proving that a completely ordered state cannot exist unless T = 0.

Another interesting point is related to the free energy difference between a mixed and a phase separated state at the same temperature and, therefore, at the same g. This free energy difference is calculated as
Δf(ϕ,̄)=f(ϕ,̄)ϕf(1,̄)(1ϕ)f(0,̄)
(67)
By using Eq. (61), we get
βΔf(ϕ,̄)=dβϵiϕ(1ϕ)+1̄ϕlnϕ+(1ϕ)ln(1ϕ),
(68)
which is essentially identical to the result of the Flory–Huggins (FH) model.6,33,47,57
Equation (68) deserves some comments. The energetic term in the original FH model6,57 also includes a (Flory) parameter accounting for the polymer-solvent and the solvent-solvent interaction. Within our field theory, it is not difficult to account too for the polymer-solvent interactions by modifying Eq. (39) as follows
B[{φμ(x)}]=κe2(d1)βϵmseβχ2×(1eβϵa)μ=1d|φμ(x)|2+eβϵaμ=1dφμ(x)2+κηe(2d1)βϵmseβχ/2μ=1dφμ1(x).
(69)
In Eq. (69), we have introduced ϵms as the interaction parameter between monomer and solvent and we have replaced ϵi with the combination χ = βϵi + 2βϵms, which closely resembles the so called Flory parameter in the original FH formulation.6,57 By repeating the exact same procedure of Sec. IV and by performing the same Legendre transform already discussed in this section, one finds
βΔf(ϕ,̄)=dχϕ(1ϕ)+1̄ϕlnϕ+(1ϕ)ln(1ϕ),
(70)
where now the energetic term includes also an explicit polymer-solvent interaction.

Beyond their formal resemblance, it must be also stressed that Eq. (70) is more general than the original FH theory, since it includes also the case when the system is polydisperse. The monodisperse limiting case is selected when dividing by the average chain length ̄ in the second term. Finally, it is evident from Eq. (70) that Δf does not depend on ϵa. Although this may appear surprising at a first sight, it is a natural consequence of the fact that all the terms featuring the bending stiffness are linear in ϕ, therefore they disappear when one computes the free energy variation as defined in Eq. (67).

The goal of the present study was to shed some light on the challenging problem of predicting the phase behavior of a system of interacting semiflexible polymers in solution because of its important consequences on protein aggregations58,59 as well as on polymer crystallization.9 

To this aim, we made the following approximations. First, we work with implicit solvent where both intrachain and interchain interactions display a short-range attraction mimicking the effect of the solvent. Second, we work on a d-dimensional lattice where these attractive energies are taken to be equal and acting only between non-consecutive nearest-neighbours chain points, and bending rigidity is represented by an energy penalty attributed on each turn of a chain. Hence the system of interacting semiflexible polymers is represented by a system of self-avoiding walks where each turn is penalized and each nearest-neighboring occurrence is rewarded. Finally, we constructed a field theory representation of this system and solved it within a mean-field approximation. Specifically, we have derived a mean-field solution of the grand potential Ω(κ, η, ϵa, ϵi) of the system [Eq. (44)] obtained by making use of a field-theoretical representation based on the polymer-magnet analogy [O(n0)-model4,13,34].

By solving the saddle point Eqs. (42) and (43), we have deduced the bond [Eq. (46)], chain [Eq. (47)] and, hence, monomer [Eq. (48)] density as a function of the parameters of the model. A discontinuity appears only (Secs. VI C and VI D) for non-zero attractive interaction between non-consecutive nearest-neighbouring monomers where there is an abrupt shift in the total monomer density ϕ (taken here as the order parameter) from a gas phase (ϕ ≃ 0) to a liquid-like phase (ϕ ≃ 1). Notably, Grand Canonical Monte Carlo simulations of the lattice model were found in very good agreement with (and, hence, confirm) the mean-field results. Notice that a similar gas-liquid transition was observed also in previous MC simulations of multi-chain systems16,17 and in experiments.9 Last but not least, our theory predicts (Sec. VII) that the free energy variation upon mixing has a Flory–Huggins-like form, yet our result is slightly more general as it accounts for the situation where the system remains polydisperse. Otherwise, since each contribution to the free energy depending on the bending stiffness is linear in ϕ, the free energy variation upon mixing is independent of the bending stiffness.

In principle, different transitions might be expected.17 First, a gas-liquid transition from a low density to a high density phase. Second, a isotropic-nematic transition from a randomly oriented isotropic phase to a phase where the stiffness induces an overall tendency to align along a common director. Finally, a coexistence of these two might also be present as the isotropic-nematic transition can be located either on the high-density (liquid) side of the gas-liquid transition for small stiffness, or on the low-density region of the low-density (gas) side for sufficiently stiff chains. Triple points then might also be present.

Quite surprisingly, the bending rigidity plays no role in our theory, its effect being to renormalize the bond and chain fugacity in apparent contradiction with numerical simulations17 predicting instead the liquid phase to be nematic. The origin of this discrepancy is likely to be ascribed to the fact that within our approach chains are polydisperse, and polydispersity is known to destabilize the isotropic-nematic transition.60 

Our approach extends previous mean-field analysis for a single chain13 to multi-chain systems by accounting, in an intuitive and transparent manner, for all the fundamental ingredients of polymer solutions (chain connectivity, bending stiffness, monomer-monomer interactions, role of dilution). While few mean-field theories do exist in the literature for multi-chain systems,33,35,37,41 none of them considered the general approach proposed here. Overall, we trust that our theory will provide a guidance toward more dedicated approaches dealing with specific cases. For instance, a recent study61 discussed the possibility of tailoring the conditions for observing the folding of a double helix from the self-folding of a single semiflexible polymer. Providing a guidance for navigating in the large parameter space that is usually required for such investigations would prove an invaluable tool.

A.G. acknowledges the MIUR PRIN-COFIN2017 Soft Adaptive Networks Grant No. 2022JWAF7Y. A.G. and A.R. acknowledge networking support by the COST Action Grant No. CA17139 (EUTOPIA).

The authors have no conflicts to disclose.

Davide Marcato: Conceptualization (equal); Formal analysis (equal); Investigation (lead); Methodology (equal); Validation (lead); Writing – original draft (equal); Writing – review & editing (equal). Achille Giacometti: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Amos Maritan: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Angelo Rosa: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In order to demonstrate Eq. (9), we start by characterizing first the consequences of the geometrical average Eq. (7) and, in particular, of the quadratic form [Eq. (8)] of the corresponding moment-generating function in the limit n → 0. For instance, it is easy to see that the m-th moment of the spin component Si,
limn0(Si)mΩ=limn0m(φi)meSφΩφ=0,
(A1)
is different from 0 only if m = 0 or m = 2. Similarly, the moment containing different components
limn0SiSjΩ=limn0φiφjeSφΩφ=0
(A2)
gives δij. To summarize, given the previous two equations it is not difficult to see that in the limit n → 0, the geometrical average Eq. (7) has the following properties:
1Ω=1,
(A3)
SiΩ=0,
(A4)
SiSjΩ=δij,
(A5)
Si1Si2SikΩ=0,ifk3.
(A6)
We focus now on Eq. (9). By using the expression for the lattice Hamiltonian Eq. (5) we have:
limn0eβHΩ=limn0xμ=1d1+JS(x)S(x+êμ)+J22(S(x)S(x+êμ))2Ω.
(A7)
It is, in fact, sufficient to stop at the second term in the expansion because of the property Eq. (A6), or every term in which a spin at a given site appears more than twice gives contribution = 0.

Interestingly, we can assign the following meaning to each of the terms in Eq. (A7):

  • The term 1 corresponds to an empty site;

  • The term S(x)S(x+êμ) corresponds to a bond connecting sites x and x+êμ;

  • The term (S(x)S(x+êμ))2 correspond to a two-step closed loop between sites x and x+êμ.

Let us consider the particular bond configuration depicted in panel (a) of Fig. 8. The corresponding term is
J4(S(x2)S(x6))(S(x6)S(x7))(S(x2)S(x3))(S(x3)S(x7))Ω=J4i1,i2,i3,i4Si1(x2)Si1(x6)Si2(x6)Si2(x7)Si3(x2)Si3(x3)Si4(x3)Si4(x7)Ω
(A8)
that, since spins on different sites are independent under the geometrical average ⟨·⟩Ω, can be factorized as
J4i1,i2,i3,i4Si1(x2)Si3(x2)ΩSi1(x6)Si2(x6)Ω×Si3(x3)Si4(x3)ΩSi2(x7)Si4(x7)Ω.
(A9)
It is now easy to realize that, in order for this term to be non-zero, the only possibility is to take i1 = i2 = i3 = i4 = i, leading to the result:
J4i=1n1=nJ4.
(A10)
It is not difficult to extend this result and conclude that every possible self-avoiding loop of k steps will appear in the expansion with a weight nJk. Let us now consider an open chain as, for instance, the one in panel (b) of Fig. 8 corresponding to the term:
J3(S(x2)S(x6))(S(x2)S(x3))(S(x3)S(x7))Ω=J3i1,i2,i3Si1(x2)Si1(x6)Si2(x2)Si2(x3)Si3(x3)Si3(x7)Ω
(A11)
Again, we can factorize
J3i1,i2,i3Si1(x2)Si2(x2)ΩSi2(x3)Si3(x3)ΩSi1(x6)ΩSi3(x7)Ω.
(A12)
Notice that, since the spins in positions x6 and x7 appear only once, because of the trace properties (A3)(A6) the weight of this configuration is = 0. We conclude then, that a single-chain configuration has non-zero weight if and only if each spin (i.e., each lattice site) appears exactly twice or does not appear at all, in other words if the configuration corresponds to a self-avoiding closed loop. Finally, let us consider the last scenario illustrated in panel (c) of Fig. 8, namely two disjointed loops. The corresponding contribution to the partition function has a similar form as of Eq. (A9), with eight occupied lattice sites instead of 4 and where again each lattice sites appears exactly twice. One can easily check that, of the eight component indices, only the indices of spins belonging to the same connected part of the configuration need to be equal in order for the term to give a non-zero contribution. In turns, this leads to the following contribution to the partition function:
J8i,j=1n1=n2J8.
(A13)
More in general, configurations with k disconnected loops have a weight proportional to nk. This concludes the proof of Eq. (9).
FIG. 8.

Examples of some possible configurations on the 2d square lattice. (a) A single closed loop. (b) An open chain. (c) A configuration made of two disjointed closed loops.

FIG. 8.

Examples of some possible configurations on the 2d square lattice. (a) A single closed loop. (b) An open chain. (c) A configuration made of two disjointed closed loops.

Close modal

Let us now see briefly why the definitions (15)(18) lead directly to the enumeration of Hamiltonian paths. The crucial point is that ⟨1⟩0 = 0. Since spins on different sites are independent under the trace operation, the only non-zero terms in Eq. (A7) are those where every spin appears exactly twice. Based on the previous discussion, it is easy to see that such terms correspond to Hamiltonian closed paths. Again, these terms will have a weight proportional to n.

To conclude, we discuss briefly the trick to count multiple open chains instead of single closed loops. The basic idea is to introduce an external magnetic field in an arbitrary direction in the O(n)-model Hamiltonian. Let us denote by H the external magnetic field along the direction 1. The expansion of eβHΩ can be truncated at
limn0x1+HS1(x)μ=1d1+JS(x)S(x+êμ)Ω.
(A14)
The term HS1(x) now corresponds to the presence of a chain end located at site x. Thus, for instance, the configuration of Fig. 8(b) would be described by the term
H2J3S1(x6)(S(x6)S(x2))(S(x2)S(x3))×(S(x3)S(x7))S1(x7)Ω.
(A15)
Notice that now, with the introduction of the external field, each spin appears exactly twice and, therefore, the weight of the configuration is non-zero. By proceeding with the factorization we get
H2J3i1,i2,i3S1(x6)Si1(x6)ΩSi1(x2)Si2(x2)Ω×Si2(x3)Si3(x3)ΩSi3(x7)S1(x7)Ω.
(A16)
The fact that the component 1 appears now explicitly determines a crucial difference with respect to the previous cases: in order for this term to be non-zero, all the component indices must be equal to 1, i.e., the direction of the external magnetic field. Therefore, only one term in the summation survives and the weight of the configuration is simply H2J3, i.e., it is of order n0. It can be easily checked with the same type of calculations that the weight of a configurations with p chains and k total bonds is H2pJk. Once again, the weight of configurations consisting of closed loops remains non-zero and of order n. Therefore, in the limit n → 0, we account for only those configurations with no contribution from closed loops.
By noticing that the denominator on the rhs of the saddle-point Eq. (49) is equal to the argument of the logarithm in the grand potential (44) and, therefore, it must be strictly positive, we may rearrange the (49) in the cubic form:
p(φ)d2κ4φ3+dκη2φ2+12dκφκη=0.
(B1)
Then we notice that the coefficient of the cubic term is positive, so p(φ → −) → − and p(φ → +) → +, while p(0) ≤ 0. In order to gain some insight on the possible solutions of (B1), we study the first derivative of p(φ) with respect to φ,
p(φ)=3d2κ4φ2+dκηφ+12dκ.
(B2)
By setting p′(φ) = 0, only two scenarios are possible (see Fig. 9):
  • If η<32 and κ<1d12η23 then p(φ) is a monotonically increasing function. Since p(0) ≤ 0, Eq. (B1) has only one real solution, >0 for any κ > 0 and = 0 for κ = 0.

  • For all other values of κ and η, p′(φ) = 0 has two solutions,
    φ±=dκη1±1+dκ123η2,
    (B3)
    corresponding to a local maximum (φ+) and a local minimum (φ) for p(φ). Since the local maximum is for φ+ < 0, and always remembering that p(0) ≤ 0, also in this case Eq. (B1) has only one real positive solution.
FIG. 9.

Two possible scenarios for the function p(φ) in Eq. (B1). In the first scenario (blue line) p(φ) is monotonically increasing, whereas in the second one (orange line) a local maximum and a local minimum appear. In both situations, the curve intersects the positive φ-semiaxis once and only once.

FIG. 9.

Two possible scenarios for the function p(φ) in Eq. (B1). In the first scenario (blue line) p(φ) is monotonically increasing, whereas in the second one (orange line) a local maximum and a local minimum appear. In both situations, the curve intersects the positive φ-semiaxis once and only once.

Close modal
It is not difficult to see that Eqs. (56) and (57) can be recast as the following two equations, both expressing ψ as a function of φ only:
ψ2dβϵi=dφ24φ+2ηdκeβϵi/2φ+ηdκeβϵi/2,
(C1)
ψ2dβϵi=12dβϵilogφ/2d2κeβϵiφ3/4dκηeβϵi/2φ2/2+dκeβϵiφ+κηeβϵi/2.
(C2)
In fact Eq. (C1) follows from the ratio of the original two, while Eq. (C2) can be obtained from Eq. (56) by extracting ψ as a function of φ. Then, it is an elementary exercise to solve numerically the system of these two equations for given values of d, η, κ and ϵi.

As an example – and without lack of generality – we consider again (Sec. VI C) the three-dimensional case study with ϵi/κBT = 1 and chain fugacity η = 0.15 (red lines in Fig. 5). Then, based on the value of the bond fugacity κ, four scenarios are possible (see Fig. 10). For low κ [panel (a)] only one solution exists, corresponding to the gas-like phase (namely, ϕbϕc ≪ 1). At some intermediate κ [panel (b)], another solution appears yet the most stable one {i.e., with the lowest grand potential βΩ [Eq. (44)]} remains the gas-like phase. Finally [panel (d)], for κ larger than some “critical” κcr [panel (c)] the most stable solution corresponds to the liquid-like one (namely, ϕ = ϕb + ϕc ≃ 1). For other chain fugacities (as well as in the case of polymers with non-zero bending stiffness, Sec. VI D) the picture remains essentially the same, thus justifying Figs. 5 and 6.

FIG. 10.

Numerical solutions (large dots) for the system of Eq. (C1) (solid lines) and Eq. (C2) (dashed lines) in the physically admissible interval φ[0,4/d] [see Eq. (46)] for d = 3, η = 0.15 and ϵi/kBT = 1 (see Sec. VI C and Fig. 5) and four representative values of the bond fugacity κ [see legends, notice the “critical” value κcr of panel (c)]. The corresponding values for the grand potential per lattice site βΩ [Eq. (44)] are also shown.

FIG. 10.

Numerical solutions (large dots) for the system of Eq. (C1) (solid lines) and Eq. (C2) (dashed lines) in the physically admissible interval φ[0,4/d] [see Eq. (46)] for d = 3, η = 0.15 and ϵi/kBT = 1 (see Sec. VI C and Fig. 5) and four representative values of the bond fugacity κ [see legends, notice the “critical” value κcr of panel (c)]. The corresponding values for the grand potential per lattice site βΩ [Eq. (44)] are also shown.

Close modal
1.
J. P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
(
Elsevier Science
,
2006
).
2.
P. N.
Pusey
and
W.
van Megen
, “
Phase behaviour of concentrated suspensions of nearly hard colloidal spheres
,”
Nature
320
(
6060
),
340
342
(
1986
).
3.
P. J.
Flory
,
Principles of Polymer Chemistry
(
Cornell University Press
,
Ithaca, NY
,
1953
).
4.
P. G.
de Gennes
,
Scaling Concepts in Polymer Physics
(
Cornell University Press
,
1979
).
5.
M.
Doi
and
S. F.
Edwards
,
The Theory of Polymer Dynamics
,
International Series of Monographs on Physics
(
Clarendon Press
,
1988
), p.
11
.
6.
M.
Rubinsteini
and
R. H.
Colby
,
Polymer Physics
(
Oxford University Press
,
Oxford
,
2003
).
7.
P.
Kawak
,
D. S.
Banks
, and
D. R.
Tree
, “
Semiflexible oligomers crystallize via a cooperative phase transition
,”
J. Chem. Phys.
155
(
21
),
214902
(
2021
).
8.
S.
Wilken
,
A.
Chaderjian
, and
O. A.
Saleh
, “
Spatial organization of phase-separated DNA droplets
,”
Phys. Rev. X
13
,
031014
(
2023
).
9.
P. D.
Olmsted
,
W. C. K.
Poon
,
T. C. B.
McLeish
,
N. J.
Terrill
, and
A. J.
Ryan
, “
Spinodal-assisted crystallization in polymer melts
,”
Phys. Rev. Lett.
81
(
2
),
373
(
1998
).
10.
J.
Bascle
,
T.
Garel
, and
H.
Orland
, “
Mean-field theory of polymer melting
,”
J. Phys. A: Math. Gen.
25
(
23
),
L1323
(
1992
).
11.
J.
Bascle
,
T.
Garel
, and
H.
Orland
, “
Some physical approaches to protein folding
,”
J. Phys. I
3
(
2
),
259
275
(
1993
).
12.
J.
Bascle
,
T.
Garel
, and
H.
Orland
, “
Formation and stability of secondary structures in globular proteins
,”
J. Phys. II
3
(
2
),
245
253
(
1993
).
13.
S.
Doniach
,
T.
Garel
, and
H.
Orland
, “
Phase diagram of a semiflexible polymer chain in a θ solvent: Application to protein folding
,”
J. Chem. Phys.
105
(
4
),
1601
1608
(
1996
).
14.
S.
Ranganathan
and
E.
Shakhnovich
, “
The physics of liquid-to-solid transitions in multi-domain protein condensates
,”
Biophys. J.
121
(
14
),
2751
2766
(
2022
).
15.
Y. J.
Sheng
,
A. Z.
Panagiotopoulos
,
S. K.
Kumar
, and
I.
Szleifer
, “
Monte Carlo calculation of phase equilibria for a bead-spring polymeric model
,”
Macromolecules
27
(
2
),
400
406
(
1994
).
16.
Y.-J.
Sheng
,
A. Z.
Panagiotopoulos
, and
S. K.
Kumar
, “
Effect of chain stiffness on polymer phase behavior
,”
Macromolecules
29
(
12
),
4444
4446
(
1996
).
17.
V. A.
Ivanov
,
M. R.
Stukan
,
M.
Muller
,
W.
Paul
, and
K.
Binder
, “
Phase diagram of solutions of stiff-chain macromolecules: A Monte Carlo simulation
,”
J. Chem. Phys.
118
(
22
),
10333
(
2003
).
18.
S. V.
Bobbili
and
S. T.
Milner
, “
Simulation study of entanglement in semiflexible polymer melts and solutions
,”
Macromolecules
53
(
10
),
3861
3872
(
2020
).
19.
R.
Everaers
,
H. A.
Karimi-Varzaneh
,
F.
Fleck
,
N.
Hojdis
, and
C.
Svaneborg
, “
Kremer–Grest models for commodity polymer melts: Linking theory, experiment, and simulation at the Kuhn scale
,”
Macromolecules
53
(
6
),
1901
1916
(
2020
).
20.
J. D.
Dietz
and
R. S.
Hoy
, “
Facile equilibration of well-entangled semiflexible bead–spring polymer melts
,”
J. Chem. Phys.
156
(
1
),
014103
(
2022
).
21.
P. I.
Kos
,
V. A.
Ivanov
, and
A. V.
Chertovich
, “
Crystallization of semiflexible polymers in melts and solutions
,”
Soft Matter
17
(
9
),
2392
2403
(
2021
).
22.
A.
Montesi
,
M.
Pasquali
, and
F. C.
MacKintosh
, “
Collapse of a semiflexible polymer in poor solvent
,”
Phys. Rev. E
69
(
2
),
021916
(
2004
).
23.
D. T.
Seaton
,
S.
Schnabel
,
D. P.
Landau
, and
M.
Bachmann
, “
From flexible to stiff: Systematic analysis of structural phases for single semiflexible polymers
,”
Phys. Rev. Lett.
110
(
2
),
028103
(
2013
).
24.
A.
Chertovich
and
P.
Kos
, “
Crumpled globule formation during collapse of a long flexible and semiflexible polymer in poor solvent
,”
J. Chem. Phys.
141
(
13
),
134903
(
2014
).
25.
M.
Marenz
and
W.
Janke
, “
Knots as a topological order parameter for semiflexible polymers
,”
Phys. Rev. Lett.
116
(
12
),
128301
(
2016
).
26.
B.
Werlich
,
M. P.
Taylor
,
T.
Shakirov
, and
W.
Paul
, “
On the pseudo phase diagram of single semi-flexible polymer chains: A flat-histogram Monte Carlo study
,”
Polymers
9
(
2
),
38
(
2017
).
27.
T.
Škrbić
,
J. R.
Banavar
, and
A.
Giacometti
, “
Chain stiffness bridges conventional polymer and bio-molecular phases
,”
J. Chem. Phys.
151
(
17
),
174901
(
2019
).
28.
T. X.
Hoang
,
A.
Giacometti
,
R.
Podgornik
,
N. T. T.
Nguyen
,
J. R.
Banavar
, and
A.
Maritan
, “
From toroidal to rod-like condensates of semiflexible polymers
,”
J. Chem. Phys.
140
(
6
),
064902
(
2014
).
29.
T. X.
Hoang
,
H. L.
Trinh
,
A.
Giacometti
,
R.
Podgornik
,
J. R.
Banavar
, and
A.
Maritan
, “
Phase diagram of the ground states of DNA condensates
,”
Phys. Rev. E
92
(
6
),
060701
(
2015
).
30.
D.
Aierken
and
M.
Bachmann
, “
Stable intermediate phase of secondary structures for semiflexible polymers
,”
Phys. Rev. E
107
(
3
),
L032501
(
2023
).
31.
U.
Bastolla
and
P.
Grassberger
, “
Phase transitions of single semistiff polymer chains
,”
J. Stat. Phys.
89
,
1061
1078
(
1997
).
32.
J. P. K.
Doye
,
R. P.
Sear
, and
D.
Frenkel
, “
The effect of chain stiffness on the phase behaviour of isolated homopolymers
,”
J. Chem. Phys.
108
(
5
),
2134
2142
(
1998
).
33.
P. J.
Flory
, “
Statistical thermodynamics of semi-flexible chain molecules
,”
Proc. R. Soc. London, Ser. A
234
(
1196
),
60
73
(
1956
).
34.
P. G.
de Gennes
, “
Exponents for the excluded volume problem as derived by the Wilson method
,”
Phys. Lett. A
38
(
5
),
339
340
(
1972
).
35.
J.
des Cloizeaux
, “
The Lagrangian theory of polymer solutions at intermediate concentrations
,”
J. Phys.
36
(
4
),
281
291
(
1975
).
36.
P. D.
Gujrati
, “
On the absence of the completely ordered phase in the flory model of semi-flexible linear polymers
,”
J. Phys. A: Math. Gen.
13
(
12
),
L437
(
1980
).
37.
P. D.
Gujrati
, “
Magnetic analog of self-avoiding random walks (polymer chains) on a lattice
,”
Phys. Rev. A
24
,
2096
2108
(
1981
).
38.
P. D.
Gujrati
and
M.
Goldstein
, “
On the validity of the Flory–Huggins approximation for semiflexible chains
,”
J. Chem. Phys.
74
(
4
),
2596
2603
(
1981
).
39.
J. C.
Wheeler
and
P.
Pfeuty
, “
The n → 0 vector model and equilibrium polymerization
,”
Phys. Rev. A
24
,
1050
1062
(
1981
).
40.
P. D.
Gujrati
, “
Lower bounds on entropy for polymer chains on a square and a cubic lattice
,”
J. Stat. Phys.
28
(
3
),
441
472
(
1982
).
41.
K. F.
Freed
, “
New lattice model for interacting, avoiding polymers with controlled length distribution
,”
J. Phys. A: Math. Gen.
18
,
871
(
1985
).
42.
H.
Orland
,
C.
Itzykson
, and
C.
de Dominicis
, “
An evaluation of the number of Hamiltonian paths
,”
J. Phys. Lett.
46
(
8
),
353
357
(
1985
).
43.
J. L.
Jacobsen
and
J.
Kondev
, “
Field theory of compact polymers on the square lattice
,”
Nucl. Phys. B
532
(
3
),
635
688
(
1998
).
44.
S.
Higuchi
, “
Field theoretic approach to the counting problem of Hamiltonian cycles of graphs
,”
Phys. Rev. E
58
,
128
132
(
1998
).
45.
A. G.
Zilman
and
S. A.
Safran
, “
Thermodynamics and structure of self-assembled networks
,”
Phys. Rev. E
66
,
051107
(
2002
).
46.
P. J.
Flory
, “
Thermodynamics of high polymer solutions
,”
J. Chem. Phys.
10
,
51
61
(
1942
).
47.
M. L.
Huggins
, “
Some properties of solutions of long-chain compounds
,”
J. Phys. Chem.
46
(
1
),
151
158
(
1942
).
48.
J.
Hubbard
, “
Calculation of partition functions
,”
Phys. Rev. Lett.
3
,
77
78
(
1959
).
49.
P. M.
Chaikin
and
T. C.
Lubensky
,
Principles of Condensed Matter Physics
(
Cambridge University Press
,
2000
).
50.
S.
Franz
,
T.
Garel
, and
H.
Orland
, “
Overlap properties and adsorption transition of two Hamiltonian paths
,”
Eur. Phys. J. B
11
(
3
),
463
468
(
1999
).
51.
S.
Lise
, “
Selected topics in polymer physics
,” Ph.D. thesis,
Scuola Internazionale Superiore di Studi Avanzati (SISSA)
,
Trieste, Italy
,
1998
.
52.
The reason to adopt “η2” for the chain fugacity is because each linear chains possesses two free ends, each of which may act as a polymerization unit.
53.
For the parsimony of notation, we adopt the same symbols “φ” and “ψ” also for the solutions of MF Eqs. (42) and (43).
54.
N.
Metropolis
,
A. W.
Rosenbluth
,
M. N.
Rosenbluth
,
A. H.
Teller
, and
E.
Teller
, “
Equation of state calculations by fast computing machines
,”
J. Chem. Phys.
21
(
6
),
1087
1092
(
1953
).
55.
M. E. J.
Newman
and
G. T.
Barkema
,
Monte Carlo Methods in Statistical Physics
(
Oxford University Press
,
Oxford
,
1999
).
56.
A preliminary analysis (data not shown) of polymer conformations generated by our GCMC algorithm seems to suggest that, on the contrary, g depends on ϕ and that, therefore, both Eqs. (65) and (66) are qualitatively wrong. At the same time, a more systematic comparison between theory and simulations would require using ϕ and ̄ as tunable parameters and, hence, a different numerical scheme. This goes beyond the present scopes and will be the subject of future research work.
57.
M.
Doi
,
Soft Matter Physics
(
Oxford University Press
,
Oxford
,
2013
).
58.
M.
Fuxreiter
and
M.
Vendruscolo
, “
Generic nature of the condensed states of proteins
,”
Nat. Cell Biol.
23
(
6
),
587
594
(
2021
).
59.
T. C. T.
Michaels
,
D.
Qian
,
A.
Šarić
,
M.
Vendruscolo
,
S.
Linse
, and
T. P. J.
Knowles
, “
Amyloid formation as a protein phase transition
,”
Nat. Rev. Phys.
5
(
7
),
379
397
(
2023
).
60.
P.
Woolston
and
J. S.
van Duijneveldt
, “
Isotropic-nematic phase transition of polydisperse clay rods
,”
J. Chem. Phys.
142
(
18
),
184901
(
2015
).
61.
J.
Du
,
H.
Yin
,
H.
Zhu
,
T.
Wan
,
B.
Wang
,
H.
Qi
,
Y.
Lu
,
L.
Dai
, and
T.
Chen
, “
Forming a double-helix phase of single polymer chains by the cooperation between local structure and nonlocal attraction
,”
Phys. Rev. Lett.
128
(
19
),
197801
(
2022
).