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–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 account^{3–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–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–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–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 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 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 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 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}

*n*-component vector,

*n*-vector model in zero external field reads

*β*= 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

*μ*”-direction. Then, we introduce the integration measure,

*δ*-function enforcing the constraint Eq. (2) and the related

*geometrical*average (or,

*trace*operation)

*n*→ 0, the following equality for the moment-generating function holds:

^{4}

*n*→ 0, the following expansion holds:

*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

^{−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 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

*Z*

_{N}in Eq. (9) – which coincides with the

*total*number of HR on our lattice – can be formally

^{50}obtained as

*Z*

_{N},

*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.

*a priori*a trace operation – that we denote by the symbol ⟨·⟩

_{0}– characterized by the desired mathematical properties:

*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).

*J*= 1 in the Hamiltonian of the

*n*-vector model [Eq. (3)], the partition function (14) of the HR is equivalent to:

*J*= 1], we have

## 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 – 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.

*grand canonical*ensemble where neither the number of chains nor the number of bonds are fixed and we introduce the grand canonical partition function,

*κ*; (b) $Nc(C)$, the total number of chains with corresponding chain fugacity

^{52}

*η*

^{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/(

*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.

*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$,

^{48,49}of the exponential of the rhs of Eq. (24) is equivalent to the expression containing the scalar field $\psi =\psi (x\u20d7)$

*d*distinct

*n*-component vectors, $S\u20d7\mu (x\u20d7)\u2261(S\mu 1(x\u20d7),S\mu 2(x\u20d7),\u2026,S\mu n(x\u20d7))$ with

*μ*= 1, 2, …,

*d*, obeying the generalized

*trace*rules:

*S*-vectors on different sites are independent from each other under the trace operations just defined.

*μ*”-direction introduced in Sec. II and

*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

*n*and they disappear in the

*n*→ 0 limit and where, by analogy with the matrix Δ [Eq. (4)], we have defined

^{48,49}of the exponential term of the last line of Eq. (32) containing the $S\u20d7$-vectors,

*d*real

*vector*fields $\phi \u20d7\mu (x\u20d7)$ (

*μ*= 1, …,

*d*) with $\phi \u20d7\mu (x\u20d7)\u2261(\phi \mu 1(x\u20d7),\phi \mu 2(x\u20d7),\u2026,\phi \mu n(x\u20d7))$ and the corresponding measure [see Eq. (26), for analogy]

*n*→ 0, one can show that, up to an unimportant multiplicative constant,

*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}

*O*(

*n*) symmetry of the vector field so that

*μ*, thus obtaining

^{13}

^{,}$q(\beta )=2+2(d\u22121)e\u2212\beta \u03f5a$. 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

*n*disappears, and thus the limit

*n*→ 0 is trivial.

## 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 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).

*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 accept

^{54}the new conformation with probability given by the expression:

*ϕ*

_{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 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,(46)$\varphi b\u2261\u27e8Nb\u27e9N=\u2212\beta \kappa \u2202\Omega (\kappa ,\eta ,\u03f5i,\u03f5a)\u2202\kappa =d4\phi 2;$ - The
*chain*density,(47)$\varphi c\u2261\u27e8Nc\u27e9N=\u2212\beta \eta 2\u2202\Omega (\kappa ,\eta ,\u03f5i,\u03f5a)\u2202\eta =d2e\beta \u03f5i\psi \kappa \eta e\u2212\beta \u03f5i/2\phi 1+de\beta \u03f5i\psi \kappa e\u2212\beta \u03f5iq(\beta )4\phi 2+\kappa \eta e\u2212\beta \u03f5i/2\phi ;$

### A. Case *ϵ*_{a} = *ϵ*_{i} = 0

*φ*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.

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 *κ* ≃ *κ**.

*ϕ*

_{b}and

*ϕ*

_{c}in the limit

*κ*≫ 1. In this regime, Eq. (49) depends only on one parameter, namely the ratio $\eta /\kappa $, and has the following solutions

*ϕ*

_{b}[Eq. (46)] and

*ϕ*

_{c}[Eq. (47)] we find

*κ*≫ 1 the monomer density

*ϕ*=

*ϕ*

_{b}+

*ϕ*

_{c}is always ≃1. However, there are two different scenarios: if $\eta /\kappa \u226a1$ 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 $\eta /\kappa \u226b1$ 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.

### B. Case *ϵ*_{a} > 0, *ϵ*_{i} = 0

*κ*→

*κq*(

*β*)/2

*d*and $\eta \u2192\eta 2dq(\beta )$. Since 2 ≤

*q*(

*β*) ≤ 2

*d*, 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

*ψ*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\u2264\psi \u22642d\beta \u03f5i$, 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*(*β*)/2*d* and $\eta \u2192\eta 2dq(\beta )$] 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} = *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. 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.

## 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.

*φ*and

*ψ*, one can express the fugacities

*η*and

*κ*in terms of

*ϕ*

_{b}and

*ϕ*

_{c}as:

*f*depends on both,

*ϕ*

_{b}and

*ϕ*

_{c}. However, it is more convenient to express it in terms of

*ϕ*and $\u2113\u0304=\varphi /\varphi c$ (i.e., a measure of the average number of monomers per chain):

*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 $\u2113\u0304\u22652$.

*U*of the system:

*not*guaranteed to be always non-negative. From Eq. (63) in fact, we readily see this to be true only when $\varphi \u22651d1\u22121\u2113\u0304$. 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\u2212(1\u22121/\u2113\u0304)$. 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 $\u2113\u0304$, we see that: (i) if $\u2113\u0304=2$ (i.e., all chains consist of 1 bond only) the total number of angles is = 0, as expected; (ii) increasing $\u2113\u0304$ 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 $\u2113\u0304$, as it should.

^{33}A central quantity of his theory is the mean bending degree of the chain,

*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 co-workers^{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.

*g*. This free energy difference is calculated as

^{6,33,47,57}

^{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

*ϵ*

_{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

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 $\u2113\u0304$ 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\u21920)$-model^{4,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 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.

## ACKNOWLEDGMENTS

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).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: DERIVATION OF EQ. (9)

*n*→ 0. For instance, it is easy to see that the

*m*-th moment of the spin component

*S*

^{i},

*m*= 0 or

*m*= 2. Similarly, the moment containing different components

*δ*

_{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:

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\u20d7(x\u20d7)\u22c5S\u20d7(x\u20d7+e\u0302\mu )$ corresponds to a bond connecting sites $x\u20d7$ and $x\u20d7+e\u0302\mu $;

The term $(S\u20d7(x\u20d7)\u22c5S\u20d7(x\u20d7+e\u0302\mu ))2$ correspond to a two-step closed loop between sites $x\u20d7$ and $x\u20d7+e\u0302\mu $.

_{Ω}, can be factorized as

*i*

_{1}=

*i*

_{2}=

*i*

_{3}=

*i*

_{4}=

*i*, leading to the result:

*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:

*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:

*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*.

*O*(

*n*)-model Hamiltonian. Let us denote by

*H*the external magnetic field along the direction 1. The expansion of $\u27e8e\u2212\beta H\u27e9\Omega $ can be truncated at

*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.

### APPENDIX B: SOLUTIONS OF EQ. (49) FOR *η* > 0

*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*′(

*φ*) = 0, only two scenarios are possible (see Fig. 9):

If $\eta <32$ and $\kappa <1d12\u2212\eta 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,corresponding to a local maximum ((B3)$\phi \xb1=\u2212d\kappa \eta 1\xb11+d\kappa \u2212123\eta 2,$*φ*_{+}) 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.

### APPENDIX C: NUMERICAL SOLUTIONS OF EQS. (56) AND (57)

*ψ*as a function of

*φ*only:

*ψ*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}/*κ*_{B}*T* = 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.

## REFERENCES

*Principles of Polymer Chemistry*

*Scaling Concepts in Polymer Physics*

*The Theory of Polymer Dynamics*

*International Series of Monographs on Physics*

*Polymer Physics*

*n*→ 0 vector model and equilibrium polymerization

*Principles of Condensed Matter Physics*

*η*

^{2}” for the chain fugacity is because each linear chains possesses two free ends, each of which may act as a polymerization unit.

*Monte Carlo Methods in Statistical Physics*

*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 $\u2113\u0304$ as tunable parameters and, hence, a different numerical scheme. This goes beyond the present scopes and will be the subject of future research work.

*Soft Matter Physics*