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 characterized1 and experimentally observed,2 in the case of polymer solutions the presence of chain connectivity3–6 makes a full understanding of its phase behavior much more challenging, in particular at high concentrations.
One of the emerging conceptual problems hinges on the difficulty to discriminate between purely kinetic effects and those associated with the underlying thermodynamics.7,8 For instance, polymers do not completely crystallize when cooled down but become structured into a hierarchy of ordered structures.9 Also, it has been argued that on cooling a polymer melt undergoes a spinodal decomposition thus making the crystallization metastable and leaving the system out of equilibrium.7 Another difficulty stems from the large number of thermodynamic and structural parameters that need to be taken into consideration: in fact, in addition to usual thermodynamic quantities such as temperature and pressure that control the system, many other microscopic parameters such as interchain (in addition to the intrachain) interactions, the number of monomers in a chain, the stiffness of the fiber and the total polymer volume fraction have to be taken into account3–6 and become axes of a large parameter space.
Particularly important appears to be the case of semiflexible polymers, as a paradigmatic example of the protein folding problem10–13 or biopolymers in a crowded environment,14 and this is currently stimulating many studies along these lines. Our current understanding of these system comes, in particular, from computer simulations which, however, have been limited so far to particular situations. A first group of studies derived the gas-liquid phase diagram for flexible15 and semiflexible16 bead-spring chains up to only 100 monomers per chain using Monte Carlo (MC) simulations. In particular, they found a phase-equilibrium diagram very similar to that of simple liquids with some minor effects ascribed to the bending rigidity. Similar results were obtained by more recent and extensive simulations.17 Other simulations aimed at understanding entanglement properties between distinct chains18–20 or the onset of nucleation process.21 However, a comprehensive picture of the phase behavior of semiflexible polymers in poor-solvent solutions is currently lacking.
Surprisingly, even in the case of a single semiflexible polymer, a general theoretical understanding of the phase behaviour is still lacking notwithstanding several studies with different techniques have recently appeared22–27 that remained, however, focused on rather specific questions. For instance, by using mean-field arguments supported by bead-spring MC simulations, it has been argued that the ground state of a single semiflexible chain can be either a rod-like or a toroidal structure depending on the bending rigidity and the contour length of the polymer28,29 and this has been confirmed recently30 by computer simulations which also accounts for the temperature dependence. Interestingly, classical studies of a lattice model31,32 observe only rod-like phases, in the form of Hamiltonian paths, likely due to the geometrical constraints imposed by the lattice. Remarkably, the observed phase diagram is in excellent agreement with mean-field predictions using a field theoretical approach.13
While a full understanding of the differences observed in lattice and off-lattice simulations is an interesting issue on its own right (and that will be discussed elsewhere), the present study will focus on providing the multi-chains counterpart of the aforementioned single chain studies.13,31,32 Specifically, in the wake of the long-standing tradition of lattice models for modeling polymer structure,10–13,33–45 here we develop a field-theoretical description of semiflexible self-avoiding chains with attractive interactions on the d-dimensional hypercubic lattice, and solve it within a mean-field approximation. Lattice grand canonical MC simulations will then be presented demonstrating the remarkable accuracy of the mean-field predictions. The multi-chain field theory approach that is proposed here build on past work by des Cloizeaux35 that extends the classical relation pointed out long ago by de Gennes4,34 of self-avoiding walks as the n → 0 limit of a spin O(n)-model, n being the number of components of each spin on a lattice point.
Notwithstanding its limitation, the present lattice approach has the great merit of making transparent the underlying physics and provide a guidance of the regions of this large parameter space that warrant a more in-depth analysis with dedicated numerical or experimental techniques. As a by-product of the theoretical analysis, we will re-obtain some classical results within a wider framework, as we will see.
The paper is organized as the following. In Sec. II, we provide a concise summary of the current understanding of the phase behavior of interacting semiflexible polymers. In particular, due to its specialized nature and for having inspired this work, we are going to highlight some salient aspects related to the past work of Orland and colleagues.10,13,42 Moreover, we will also discuss some conclusions from recent computational work by other groups which has the advantage of providing a broader view beyond the mean-field approach. The novel part of this work starts in Sec. III A, where we will introduce the exact grand canonical partition function Z of the lattice model for a multi-chain system (i.e., a polymer solution) on the d-dimensional cubic lattice which takes into account the local bending stiffness of the polymer fiber, excluded-volume and short-range attractive interactions between close-by monomers. Then, by exploiting the analogy between self-avoiding polymers and the spin O(n → 0)-model, we construct (Sec. III B) the exact field-theoretic representation of Z. Since an exact computation of Z is unfeasible, we describe a uniform saddle-point approximation (Sec. IV) and obtain the corresponding mean-field solution of the problem, the reliability of which is demonstrated by comparison to MC computer simulations (Secs. V and VI). Finally we show that our results (Sec. VII) recapitulate, as particular cases, several models that have been discussed in the past and we demonstrate the equivalence between our approach and the classical Flory–Huggins theory46,47 for mixtures. Discussion and conclusions, with an outline on open problems and possible future perspective, are presented in Sec. VIII.
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/ad = (L/a)d where V is the volume of the lattice: that is, we consider all lattice points occupied and no vacancies. For simplicity, we further restrict our considerations to closed paths, or Hamiltonian rings (HR), knowing that the statistics is the same in the thermodynamic limit.42
We demonstrate now that there exists an alternative method51 of finding ZN (14) that represents both a shortcut with respect to the approach presented so far and has the great advantage of being exportable to the more general situation considered here (Sec. III) in a relatively straightforward manner.
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 = Ld and N = V/ad. Chains are self- and mutually-avoiding, i.e., any two monomers – be they from the same or different chains – cannot occupy the same lattice site. Chain stiffness is modeled by introducing a bending energy penalty ϵa > 0 for two consecutive bonds along the same chain forming a turn (or, an angle), while attractive interactions between non-bonded monomers are accounted for by an energy reward −ϵi < 0 (ϵa, ϵi > 0) for any two monomers which are separated by a unit lattice distance and either are non-consecutive if they belong to the same chain or they are on distinct chains.
Schematic illustration of a particular configuration on the square lattice (d = 2). By assuming periodic boundary conditions (PBC, see text for details), we have: chains, bonds, turns or angles (marked as blue corners), and pairs of interacting monomers (connected by dashed red lines, with two pairs interacting through PBC.
Schematic illustration of a particular configuration on the square lattice (d = 2). By assuming periodic boundary conditions (PBC, see text for details), we have: chains, bonds, turns or angles (marked as blue corners), and pairs of interacting monomers (connected by dashed red lines, with two pairs interacting through PBC.
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.
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
V. MONTE CARLO SIMULATIONS
In order to check the validity of the MF approximation as well as to assess its limits, we have performed Metropolis54 grand canonical Monte Carlo (GCMC) computer simulations of the lattice model (Sec. III) on the three-dimensional cubic lattice. Essentially, the goal of the GCMC simulations is to obtain a representative sample of polymer configurations in agreement with the grand canonical partition function (21).
In order to check for finite-size effects, we have performed preliminary calculations and compared corresponding results for lattice sizes L/a = 4, 8, 16. This analysis indicates that for L/a ≥ 8 the results do not change significantly with the lattice size and hence we will fix L/a = 8 in all calculations henceforth. This guarantees a good compromise between computational efficiency and accuracy.
For each chosen pair of (κ, η) values, we have led a simulation run consisting of 107 MC steps. Then, for each GCMC trajectory, a standard block analysis procedure55 has been carried out in order to estimate uncertainties on the considered physical observables. Every trajectory has been checked individually in order to make sure that all the curves obtained from corresponding block analyses and representing the MC-time evolution of the distinct quantities have completely equilibrated. This procedure has been applied to the number of bonds, the number of chains and the internal energy of the system.
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)
- The chain density,(47)
A. Case ϵa = ϵi = 0
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.
ϵa = ϵi = 0. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are, respectively, for the MF solution and the GCMC computer simulations.
ϵa = ϵi = 0. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are, respectively, for the MF solution and the GCMC computer simulations.
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. Case ϵa > 0, ϵi = 0
By fixing the bending stiffness to the convenient value ϵa = kBT, MF calculations for the bond and chain density, ϕb and ϕc, as a function of the bond fugacity κ and for chain fugacities η = 0.2 and η = 1.5 are shown as solid lines in Fig. 3 [panels (a) and (b), respectively] and compared to corresponding GCMC simulations (symbols). As in previous Sec. VI A, the agreement between MF calculation and GCMC simulations is remarkable.
ϵa/kBT = 1, ϵi = 0. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are as in Fig. 2.
ϵa/kBT = 1, ϵi = 0. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are as in Fig. 2.
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).
ϵa/kBT = 1, ϵi = 0. Angle density ϕa as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are as in Fig. 2.
ϵa/kBT = 1, ϵi = 0. Angle density ϕa as a function of the bond fugacity κ and for chain fugacities η = 0.2 (blue) and η = 1.5 (red). Solid lines and symbols are as in Fig. 2.
C. Case ϵa = 0, ϵi > 0
In fact, by fixing the monomer-monomer interaction to the convenient value ϵi = kBT, MF calculations for the bond and chain density, ϕb and ϕc, as a function of the bond fugacity κ and for chain fugacities η = 0.15 and η = 0.25 are shown as solid lines in Fig. 5 [panels (a) and (b), respectively] and compared to corresponding GCMC simulations (symbols). In both panels one can easily distinguish two different phases, with both densities acting as order parameters. In the first (gas-like) phase, the total monomer density ϕ = ϕb + ϕc is close to 0 and ϕb/ϕc ≃ 1. This is valid up to some critical value κcr above which ϕ is close to 1 and ϕb/ϕc > 1 (liquid-like phase).
ϵa = 0, ϵi/kBT = 1. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.15 (blue) and η = 0.25 (red). The discontinuity predicted by MF calculations (lines) is confirmed by GCMC simulations (symbols), yet the critical value of κ is slightly different between theory and simulations.
ϵa = 0, ϵi/kBT = 1. Bond density ϕb (a) and chain density ϕc (b) as a function of the bond fugacity κ and for chain fugacities η = 0.15 (blue) and η = 0.25 (red). The discontinuity predicted by MF calculations (lines) is confirmed by GCMC simulations (symbols), yet the critical value of κ is slightly different between theory and simulations.
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 ] to absorb the terms accounting for the bending stiffness and obtain a system of equations which is equivalent to that presented in Sec. VI C. The general behavior of the solutions will thus be exactly the same (and, so, the discontinuities in ϕb and ϕc) as for a system with no bending stiffness, the only difference being in the changing of the coexistence line between the two phases in the (κ, η)-plane, depending on the value of ϵa.
By fixing again the bending stiffness and the monomer-monomer interaction to the values ϵa = ϵi = kBT, MF calculations for the bond and chain density, ϕb and ϕc, as a function of the bond fugacity κ and for chain fugacities η = 0.15 and η = 0.25 are shown as solid lines in Fig. 6 [panels (a) and (b), respectively] and compared to corresponding GCMC simulations (symbols). As expected, GCMC simulations confirm MF calculations and we can distinguish once again the gas (ϕ = ϕb + ϕc ≈ 0) and liquid (ϕ ≈ 1) phases. Finally, analogously to Sec. VI C, we produce examples of coexistence and gas-liquid transition lines, see panel (a) [dashed lines, different colors are for different ϵi (see caption)] and panel (b) [empty symbols, different colors are for different η (see caption)] of Fig. 7.
ϵ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.
ϵ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.
(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.
(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.
Note that a comparison between the two estimates gives g < gF. On the other hand we can also argue that gF is clearly overestimating the true value because it means that at each position along the chain, the possible directions available to make a turn are always 2(d − 1), thereby not accounting for the long-range correlations due to self-avoidance. By contrast, Eq. (66) implies that the number of possible directions to make a turn at each step along the chain is on average d − 1, thus effectively accounting for possible long-range correlations induced by self-avoidance.
As both estimates Eqs. (65) and (66) are independent of the concentration as remarked, it seems not plausible that it may act as an order parameter for a phase transition. Nevertheless, Flory postulates the existence of a phase-separated, ordered state (g = 0) also at T > 0, and assumes that the entropy of such state is 0. He then proceeds in computing the free energy difference between the latter and a completely mixed, disordered state (g ≠ 0), the free energy of which is computed within the same mean-field theory. He then verifies a posteriori that there exists a critical temperature below which the ordered state is thermodinamically more stable. The premise that the entropy of the ordered state is 0 is crucial for the derivation of Flory, and it is ultimately this assumption that leads to the appearance of a critical temperature, corresponding to a first-order phase transition. Later studies by Gujrati and co-workers36,38,40 were however able to derive an exact lower bound for the entropy of a system of self-avoiding chains on a lattice that was found to be strictly positive at any temperature T > 0, therefore proving that a completely ordered state cannot exist unless T = 0.
Beyond their formal resemblance, it must be also stressed that Eq. (70) is more general than the original FH theory, since it includes also the case when the system is polydisperse. The monodisperse limiting case is selected when dividing by the average chain length in the second term. Finally, it is evident from Eq. (70) that Δf does not depend on ϵa. Although this may appear surprising at a first sight, it is a natural consequence of the fact that all the terms featuring the bending stiffness are linear in ϕ, therefore they disappear when one computes the free energy variation as defined in Eq. (67).
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 aggregations58,59 as well as on polymer crystallization.9
To this aim, we made the following approximations. First, we work with implicit solvent where both intrachain and interchain interactions display a short-range attraction mimicking the effect of the solvent. Second, we work on a d-dimensional lattice where these attractive energies are taken to be equal and acting only between non-consecutive nearest-neighbours chain points, and bending rigidity is represented by an energy penalty attributed on each turn of a chain. Hence the system of interacting semiflexible polymers is represented by a system of self-avoiding walks where each turn is penalized and each nearest-neighboring occurrence is rewarded. Finally, we constructed a field theory representation of this system and solved it within a mean-field approximation. Specifically, we have derived a mean-field solution of the grand potential Ω(κ, η, ϵa, ϵi) of the system [Eq. (44)] obtained by making use of a field-theoretical representation based on the polymer-magnet analogy [-model4,13,34].
By solving the saddle point Eqs. (42) and (43), we have deduced the bond [Eq. (46)], chain [Eq. (47)] and, hence, monomer [Eq. (48)] density as a function of the parameters of the model. A discontinuity appears only (Secs. VI C and VI D) for non-zero attractive interaction between non-consecutive nearest-neighbouring monomers where there is an abrupt shift in the total monomer density ϕ (taken here as the order parameter) from a gas phase (ϕ ≃ 0) to a liquid-like phase (ϕ ≃ 1). Notably, Grand Canonical Monte Carlo simulations of the lattice model were found in very good agreement with (and, hence, confirm) the mean-field results. Notice that a similar gas-liquid transition was observed also in previous MC simulations of multi-chain systems16,17 and in experiments.9 Last but not least, our theory predicts (Sec. VII) that the free energy variation upon mixing has a Flory–Huggins-like form, yet our result is slightly more general as it accounts for the situation where the system remains polydisperse. Otherwise, since each contribution to the free energy depending on the bending stiffness is linear in ϕ, the free energy variation upon mixing is independent of the bending stiffness.
In principle, different transitions might be expected.17 First, a gas-liquid transition from a low density to a high density phase. Second, a isotropic-nematic transition from a randomly oriented isotropic phase to a phase where the stiffness induces an overall tendency to align along a common director. Finally, a coexistence of these two might also be present as the isotropic-nematic transition can be located either on the high-density (liquid) side of the gas-liquid transition for small stiffness, or on the low-density region of the low-density (gas) side for sufficiently stiff chains. Triple points then might also be present.
Quite surprisingly, the bending rigidity plays no role in our theory, its effect being to renormalize the bond and chain fugacity in apparent contradiction with numerical simulations17 predicting instead the liquid phase to be nematic. The origin of this discrepancy is likely to be ascribed to the fact that within our approach chains are polydisperse, and polydispersity is known to destabilize the isotropic-nematic transition.60
Our approach extends previous mean-field analysis for a single chain13 to multi-chain systems by accounting, in an intuitive and transparent manner, for all the fundamental ingredients of polymer solutions (chain connectivity, bending stiffness, monomer-monomer interactions, role of dilution). While few mean-field theories do exist in the literature for multi-chain systems,33,35,37,41 none of them considered the general approach proposed here. Overall, we trust that our theory will provide a guidance toward more dedicated approaches dealing with specific cases. For instance, a recent study61 discussed the possibility of tailoring the conditions for observing the folding of a double helix from the self-folding of a single semiflexible polymer. Providing a guidance for navigating in the large parameter space that is usually required for such investigations would prove an invaluable tool.
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)
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 corresponds to a bond connecting sites and ;
The term correspond to a two-step closed loop between sites and .
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.
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.
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.
APPENDIX B: SOLUTIONS OF EQ. (49) FOR η > 0
If and then p(φ) is a monotonically increasing function. Since p(0) ≤ 0, Eq. (B1) has only one real solution, for any κ > 0 and = 0 for κ = 0.
- For all other values of κ and η, p′(φ) = 0 has two solutions,corresponding to a local maximum (φ+) and a local minimum (φ−) for p(φ). Since the local maximum is for φ+ < 0, and always remembering that p(0) ≤ 0, also in this case Eq. (B1) has only one real positive solution.(B3)
Two possible scenarios for the function p(φ) in Eq. (B1). In the first scenario (blue line) p(φ) is monotonically increasing, whereas in the second one (orange line) a local maximum and a local minimum appear. In both situations, the curve intersects the positive φ-semiaxis once and only once.
Two possible scenarios for the function p(φ) in Eq. (B1). In the first scenario (blue line) p(φ) is monotonically increasing, whereas in the second one (orange line) a local maximum and a local minimum appear. In both situations, the curve intersects the positive φ-semiaxis once and only once.
APPENDIX C: NUMERICAL SOLUTIONS OF EQS. (56) AND (57)
As an example – and without lack of generality – we consider again (Sec. VI C) the three-dimensional case study with ϵi/κBT = 1 and chain fugacity η = 0.15 (red lines in Fig. 5). Then, based on the value of the bond fugacity κ, four scenarios are possible (see Fig. 10). For low κ [panel (a)] only one solution exists, corresponding to the gas-like phase (namely, ϕb ≃ ϕc ≪ 1). At some intermediate κ [panel (b)], another solution appears yet the most stable one {i.e., with the lowest grand potential βΩ [Eq. (44)]} remains the gas-like phase. Finally [panel (d)], for κ larger than some “critical” κcr [panel (c)] the most stable solution corresponds to the liquid-like one (namely, ϕ = ϕb + ϕc ≃ 1). For other chain fugacities (as well as in the case of polymers with non-zero bending stiffness, Sec. VI D) the picture remains essentially the same, thus justifying Figs. 5 and 6.
Numerical solutions (large dots) for the system of Eq. (C1) (solid lines) and Eq. (C2) (dashed lines) in the physically admissible interval [see Eq. (46)] for d = 3, η = 0.15 and ϵi/kBT = 1 (see Sec. VI C and Fig. 5) and four representative values of the bond fugacity κ [see legends, notice the “critical” value κcr of panel (c)]. The corresponding values for the grand potential per lattice site βΩ [Eq. (44)] are also shown.
Numerical solutions (large dots) for the system of Eq. (C1) (solid lines) and Eq. (C2) (dashed lines) in the physically admissible interval [see Eq. (46)] for d = 3, η = 0.15 and ϵi/kBT = 1 (see Sec. VI C and Fig. 5) and four representative values of the bond fugacity κ [see legends, notice the “critical” value κcr of panel (c)]. The corresponding values for the grand potential per lattice site βΩ [Eq. (44)] are also shown.