Phase behaviour of semiflexible lattice polymers in poor-solvent solution: mean-field theory and Monte Carlo simulations

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.


I. INTRODUCTION
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 characterized 1 and experimentally observed 2 , in the case of polymer solutions the presence of chain connectivity [3][4][5][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 a) Electronic mail: dmarcato@sissa.itb) Electronic mail: achille.giacometti@unive.itc) Electronic mail: amos.maritan@pd.infn.itd) Electronic mail: anrosa@sissa.it of the fiber and the total polymer volume fraction have to be taken into account [3][4][5][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 problem [10][11][12][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 flexible 15 and semiflexible 16 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 chains [18][19][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 appeared [22][23][24][25][26][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 polymer 28,29 and this has been confirmed recently 30 by computer simulations which also accounts for the temperature dependence.Interestingly, classical studies of a lattice model 31,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 longstanding tradition of lattice models for modeling polymer structure [10][11][12][13][33][34][35][36][37][38][39][40][41][42][43][44][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. Lattie 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 Cloizeaux 35 that extends the classical relation pointed out long ago by de Gennes 4,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 indepth 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 Section 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 theory 46,47 for mixtures.Discussion and conclusions, with an outline on open problems and possible future perspective, are presented in Sec.VIII.

II. REVIEW OF THE SINGLE-CHAIN FORMALISM: HAMILTONIAN RINGS
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 /a d = (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, is associated to each lattice point ⃗ x.By defining the scalar product ⃗ S(⃗ x) ) between any two vectors associated to lattice points ⃗ x and ⃗ x ′ , we assume the following constraint on the norm-square of ⃗ S(⃗ x): The reduced Hamiltonian for the n-vector model in zero external field reads where β = 1/(k B T ), k B being the Boltzmann constant and T being the temperature, J is the coupling constant between the spins, and where From now on we assume periodic boundary conditions and Eq. ( 3) can be also written as where êµ is the unit lattice vector pointing towards the "+µ"-direction.Then, we introduce the integration measure, with the Dirac δ-function enforcing the constraint Eq. ( 2) and the related geometrical average (or, trace operation) 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 : 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: 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 where 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, illdefined.However, this technical difficulty can be overcome through the more rigorous approach 42,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 We now notice that Z N in Eq. ( 9) -which coincides with the total number of HR on our lattice -can be formally 50 obtained as Finally, by combining Eq. ( 12) with ( 13) we get the following expression for Z N , which was introduced first in 42 .It is important to stress that, in order to compute Z N , 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 Z N , 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 method 51 of finding Z N (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: 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: In fact, by taking the Hubbard-Stratonovich transformation of the term inside brackets (Eq.(10) with J = 1), we have 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.

III. THE MANY-CHAIN FIELD THEORY A. The model
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 = L d and N = V /a d .Chains are self-and mutually-avoiding, i.e. any two monomers -be they from the same or different chains -can not 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, where the sum is intended over the set of all possible configurations {C} and where the thermal (Boltzmann) weight for each conformation C,

PBC
depends on: (a) N b (C), the total number of bonded monomer pairs with corresponding bond fugacity κ; (b) N c (C), the total number of chains with corresponding chain fugacity 52 η 2 ; (c) N a (C) and N i (C), respectively the total number of corners and the total number of nonbonded monomer pairs separated by one lattice distance.Again, β = 1/(k B T ) is the Boltzmann factor at temperature T and k B 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.

B. Partition function and field-theoretic representation
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 (N b ), chains (N c ), angles (N a ), and non-bonded monomer-monomer pairs (N i ) characterizing each given chain configuration C. To this purpose, we start by defining the scalar function for the configuration C, By using the matrix ∆(⃗ x, ⃗ x ′ ) (Eq. ( 4)), we have and the Hubbard-Stratonovich transformation 48,49 of the exponential of the r.h.s of Eq. ( 24) is equivalent to the expression containing the scalar field ψ = ψ(⃗ x) Then, at each lattice position ⃗ x we introduce d n- d, obeying the generalized trace rules: 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, where êµ is the unit lattice vector pointing towards the "+µ"-direction introduced in Sec.II and 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 n k 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 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 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 transformation 48,49 of the exponential term of the last line of Eq. ( 32) containing the ⃗ S-vectors, e where we have introduced the ) and the corresponding measure (see Eq. ( 26), for anal-ogy) 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, where we have defined the following functionals: 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.

IV. MEAN-FIELD SOLUTION: SADDLE-POINT APPROXIMATION
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 for every ⃗ x and every µ, thus obtaining where 13 q(β) = 2 + 2(d − 1)e −βϵa .In terms of the solutions 53 φ = φ(κ, η, ϵ i , ϵ a ) and ψ = ψ(κ, η, ϵ i , ϵ a ) of the MF Eqs. ( 42) and ( 43), the grand potential per lattice site (up to unimportant additive constants) reads 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 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.

V. MONTE CARLO SIMULATIONS
In order to check the validity of the MF approximation as well as to assess its limits, we have performed Metropolis 54 Grand Canonical Monte Carlo (GCMC) computer simulations of the lattice model (Sec.III) on the threedimensional 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 (N b (C)), distinct chains (N c (C)), turns (N a (C)) and pairs of non-bonded nearestneighbor monomers (N i (C)).Therefore, at each MC step one single bond is randomly inserted in or removed from the lattice, provoking a change of the configuration C 0 to the configuration C 1 .In order to enforce the condition of detailed balance, we accept 54 the new conformation with probability given by the expression: w(κ,η,ϵa,ϵi; C1) w(κ,η,ϵa,ϵi; C0) where ϕ b,0 (respectively, ϕ b,1 ) is the bond density (see Eq. ( 46)) of the configuration C 0 (resp., C 1 ).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 10 7 MC steps.Then, for each GCMC trajectory, a standard block analysis procedure 55 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.

VI. MEAN-FIELD SOLUTION VS. MONTE CARLO SIMULATIONS
From our mean-field estimate of the grand potential, Eq. ( 44), we can compute: • The bond density, • The chain density, • The total monomer density, In the rest of this Section, we specialize the saddle-point equations ( 42) and (43) to various particular cases and compare the corresponding results to Monte Carlo simulations in d = 3.
A. Case ϵa = ϵi = 0 The simplest case to be considered is the case of noninteracting (still non-overlapping) flexible chains with no bending penalty nor monomer-monomer attractive interactions.In spite its simplicity, this case proves to be rather instructive.In this case Eqs. ( 42) and (43) read The only relevant field is thus φ and, since it satisfies the simple cubic equation (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.
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 By plugging these results into the expressions for ϕ b (Eq.( 46)) and ϕ c (Eq. ( 47)) we find and 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 It is easy to see that this is the same situation of Sec.VI A with renormalized fugacities κ → κq(β)/2d and η → η 2d q(β) .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 = k B T , 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.
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).C. Case ϵa = 0, ϵi > 0 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 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 = k B T , 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).
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.D. Case ϵa > 0, ϵi > 0 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 η → η 2d q(β) ) to absorb the terms accounting for the bending stiffness and obtain a  (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.
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 monomermonomer interaction to the values ϵ a = ϵ i = k B T , 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.(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.

VII. CONNECTIONS TO THE FLORY-HUGGINS THEORY OF MIXING
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: This allows us to compute the reduced free energy per site where it is clear that f depends on both, ϕ b and ϕ c .However, it is more convenient to express it in terms of ϕ and l = ϕ/ϕ c (i.e., a measure of the average number of monomers per chain): where e = 2.71828... is the Euler's number.Eq. ( 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 l ≥ 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: where Eqs. ( 63) and (64) bear interesting physical interpretations.Let us first discuss Eq. ( 63).Within our meanfield approach the number of interactions per lattice site ⟨Ni⟩ N is not guaranteed to be always non-negative.From Eq. (63) in fact, we readily see this to be true only when 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−(1−1/ l).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 l, we see that: (i) if l = 2 (i.e., all chains consist of 1 bond only) the total number of angles is = 0, as expected; (ii) increasing l 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 l, 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, 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.
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 < g F .On the other hand we can also argue that g F 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 coworkers 36,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 By using Eq. ( 61), we get which is essentially identical to the result of the Flory-Huggins (FH) model 6,33,47,57 .
Eq. ( 68) deserves some comments.The energetic term in the original FH model 6,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 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 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 l 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).

VIII. DISCUSSION AND CONCLUSIONS
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 aggregations 58,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(n → 0)model 4,13,34 ).
By solving the saddle point equations ( 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 systems 16,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 simulations 17 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 chain 13 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 study 61 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.
It is now easy to realize that, in order for this term to be non-zero, the only possibility is to take i 1 = i 2 = i 3 = i 4 = i, leading to the result: 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 nJ k .Let us now consider an open chain as, for instance, the one in panel (b) of Fig. 8 corresponding to the term: ⟨S i1 (⃗ x 2 )S i1 (⃗ x 6 ) S i2 (⃗ x 2 )S i2 (⃗ x 3 ) S i3 (⃗ x 3 )S i3 (⃗ x 7 )⟩ Ω (A11) Again, we can factorize ⟨S i1 (⃗ x 2 )S i2 (⃗ x 2 )⟩ Ω ⟨S i2 (⃗ x 3 )S i3 (⃗ x 3 )⟩ Ω ⟨S i1 (⃗ x 6 )⟩ Ω ⟨S i3 (⃗ x 7 )⟩ Ω . (A12) Notice that, since the spins in positions ⃗ x 6 and ⃗ x 7 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 8 occupied lattice sites instead of 4 and where again each lattice sites appears exactly twice.One can easily check that, of the 8 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: More in general, configurations with k disconnected loops have a weight proportional to n k .This concludes the proof of Eq. ( 9).
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 lim (A14) The term HS 1 (⃗ 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 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 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 H 2 J 3 , i.e. it is of order n 0 .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 H 2p J k .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.44)) are also shown.

PBCFIG. 1 .
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, N b (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. 6 .
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.
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. 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. 8 .
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.