Role of interaction anisotropy in polymer cononsolvency: insights from the Flory–Huggins–Potts framework

Abstract

Cononsolvency occurs when mixing two good solvents creates poor-solvent conditions for polymers over specific composition ranges, causing macroscopic phase separation or microscopic chain collapse. Despite its technological and biophysical relevance, the connection between macroscopic and microscopic manifestations of cononsolvency remains unclear. A key challenge is identifying which interactions govern cononsolvency: coarse-grained analyses like standard Flory–Huggins models assume purely isotropic interactions, while atomistic simulations contain complex anisotropic interactions that cannot be precisely controlled or isolated. Here, we address the role of interaction anisotropy using the Flory–Huggins–Potts framework, which yields χ as a thermodynamic average over both configurational and internal-state coarse-grained degrees of freedom. This enables controlled comparison between systems with isotropic versus orientation-dependent interactions that share identical effective χ parameters, either driving cononsolvency by strong solvent–cosolvent affinity or preferential polymer–cosolvent affinity. While pairs of systems exhibit equivalent macroscopic phase behavior, lattice Monte Carlo simulations reveal that those featuring anisotropic or orientation-dependent interactions generate distinct collapse signatures, particularly in reentrant coil–globule transitions or characteristics of the solvation structure. These results demonstrate how microscopic interactions influences cononsolvency behavior beyond what effective χ parameters alone predict.

---

1 Introduction

Cononsolvency in polymer solutions arises when mixing two individually good solvents yields poor-solvent conditions for a polymer. Under such conditions, at a macroscopic level, the solution may phase separate^1–5 or, at a microscopic level, the characteristic size of polymer chains may decrease.^6–17 This composition-sensitive behavior has relevance across scientific and technological applications, such as detecting volatile organic compounds by measuring the response of a hydrogel,¹⁸ measuring enantiomeric excess in organic synthesis by observing the fluorescent enhancement,¹⁹ altering rates of organic reactions at liquid–liquid interfaces,²⁰ tuning the frictional properties of polymer brushes by adjusting solution composition,^21,22 and affecting the transport of molecules through a porous polymer matrix.^21,23,24 Furthermore, the physics governing cononsolvency has profound implications for understanding biological phase separation, protein denaturation,^25–30 and cell recovery in tissues.^31–33 There is thus significant interest to understand polymer cononsolvency and related phenomena.

Multiple mechanisms of polymer cononsolvency have been proposed on the basis of theoretical, computational, and experimental analyses. The preferential mixing mechanism attributes polymer collapse to favorable solvent–cosolvent interactions, where the polymer sacrifices conformational entropy to promote solvent–cosolvent mixing.^{4,13,16,34–36} The seminal experimental study by Schild et al.³⁷ demonstrated, however, that perturbing water–methanol interactions (χ₁₂) is insufficient to predict cononsolvency of PNIPAM (poly(N-isopropylacrylamide)) in a methanol–water solution. An alternative mechanism involves preferential polymer–cosolvent adsorption, where trace amounts of cosolvent selectively solvate the polymer and form “enthalpic bridges” between distal chain segments.³⁸ These interactions can trigger phase separation or single-chain collapse.^14,39–46 Effects akin to preferential adsorption can also arise from entropic factors, such as size asymmetry between solvent and cosolvent leading to effective depletion interactions.³⁶ In another scenario, the cosolvent can act like a surfactant, such that polymer contraction allows the cosolvent to interact favorably with both solvent and polymer.⁴⁷ The concept of geometric frustration has also been used to explain cononsolvency, proposing that competition between cosolvent and solvent for polymer solvation destabilizes the local environment.^48,49 Overall, these studies highlight polymer cononsolvency as a phenomenon with various competing physicochemical interactions.

Previous studies using Flory–Huggins (FH) theory have demonstrated that phase separation can be induced when either solvent–cosolvent or polymer–solvent interactions dominate, even when the polymer is miscible with each solvent individually.^34,50,51 Furthermore, analyses of structure factors obtained following application of the random phase approximation provide similar conclusions.⁵² Nevertheless, mean-field insights into the energetic factors driving phase separation may not always translate into understanding microscopic physics, such as single-chain conformational behavior.^38,53 In this direction, field-theoretic methods^14,46 have demonstrated that single-chain collapse indeed occurs given strong, preferential affinity of the polymer for one of the solvents, but there is no apparent coil–globule transition driven purely by solvent–cosolvent affinity. Recent developments in variational field-theoretic treatments⁵⁴ may enable further insights into single-chain behavior under such conditions.

Molecular dynamics (MD) simulations of systems with generic chemical attributes (e.g., bead-spring polymers in monomeric solvents) have effectively illustrated many of the aforementioned mechanisms of cononsolvency at a microscopic level.^36,38,55 A typical observation from such simulations is that the R_g for a single polymer chain first decreases and then increases as cosolvent is progressively added. This scenario has been induced by solvent–cosolvent mixing, preferential polymer–solvent adsorption, and the surfactant-like action of the cosolvent by analyzing the molecular environment of the polymer and carefully evaluating molecular affinities.^{38,49,56–58}

MD simulations have also been used to examine the behavior of chemically specific systems, such as PNIPAM in water and ethanol. In such cases, single-chain collapse, consistent with cononsolvency, has been observed. Despite significant interest and detailed study, there is debate as to what precise molecular interactions or combination thereof begets this behavior.^59–63 This highlights how microscopic drivers of cononsolvency can be non-trivial to resolve, even if the basic physics are well understood.

Several recent works have established quantitative connections between microscopic chain behavior and macroscopic phase separation in cononsolvency for systems with isotropic interactions using lattice models. The foundational justification for lattice models in polymer physics originates from de Gennes, who demonstrated that self-avoiding chains correspond to the n → 0 limit of -vector spin models, thereby linking polymer statistics to critical phenomena and establishing lattice theories as a rigorous framework for polymer scaling and mean-field behavior.^64–67 Zhang has comprehensively mapped the energetic conditions (i.e., χ parameters) that yield cononsolvency through mean-field analysis.^51,68 Previously, Zhang et al. also demonstrated via Wang–Landau simulations that trace cosolvent addition induces continuous coil–globule–coil transitions without requiring explicit solvent–solvent attraction.⁶⁹ By a complementary and distinct approach, Marcato et al. mapped lattice models onto -vector spin models to derive exact partition functions and field-theoretic descriptions amenable to analytics.⁷⁰ Li et al. have established connections between single-chain collapse, multi-chain aggregation, and mean-field theory for block copolymers and homopolymer chain(s) in binary solvents, in the context of FH models.⁷¹ While these studies establish that specific combinations of scalar χ parameters derived from isotropic interactions produce cononsolvency, this leaves unresolved whether orientation-dependent interactions fundamentally alters cononsolvency mechanisms or merely modulates existing ones.

While simple FH theory captures cononsolvency phenomena, the apparent connection between thermoresponsive polymer solutions and cononsolvency^2,41,72 and highlighted importance of hydrogen bonding,^41,73 ultimately suggest that orientation-dependent interactions may be relevant for at least some specific cononsolvency mechanisms. Here, we directly investigate how anisotropic interactions influences polymer cononsolvency in the context of a Flory–Huggins–Potts (FHP) framework,⁷⁴ which previously demonstrated that including orientation-dependent interactions within FH-like models enabled description of thermoresponsive phenomena (e.g., miscibility loops and heating-induced coil–globule transitions) without temperature-dependent χ parameters. This work does not resolve mechanisms or debate for any specific system; instead, it clarifies how orientation-dependent interactions, like hydrogen bonding or even packing effects, influence cononsolvency phenomena within a single, tractable framework. In particular, we demonstrate how systems with the same effective χ interaction parameters, and thus identical macroscopic phase behavior, may exhibit different microscopic signatures of cononsolvency.

The remainder of the paper is organized as follows. After presentation of methods, we first identify established phase-separation regimes driven by solvent–cosolvent mixing versus polymer–cosolvent adsorption. Then, we compare corresponding coil–globule–coil transitions across systems with identical effective χ parameters but different underlying energetic contributions using lattice Monte Carlo simulations. Subsequent analysis reveals how anisotropic interactions generate microscopically distinct collapse pathways compared to isotropic systems, despite yielding identical macroscopic χ parameters. These findings provide additional insights into cononsolvency, by ascertaining the influence of additional degrees of freedom imparted by the FHP framework, and point towards strategies for tuning polymer response based on molecular-level interactions.

2 Methods

2.1 Flory–Huggins–Potts framework and stability analysis

To enable description of orientation-dependent interactions, we extend the recently developed FHP framework⁷⁴ to model ternary polymer solutions. In effect, the FHP framework augments conventional FH theory by assigning internal states to particles that modulate pairwise interactions. These states may be interpreted, for example, as capturing aspects of molecular orientation at a coarse-grained resolution, enabling expression of orientation-dependent interactions, whereas traditional FH reflects mean-field, isotropic interactions. While the motivation and detailed development of the FHP framework can be found in ref. 74, we briefly review some salient aspects.

The behavior of the FHP model is governed by a well-defined Hamiltonian. We consider a fully occupied lattice of polymer, solvent, and cosolvent sites, with the system energy given by bonded and non-bonded contributions:

where α_i is the species type (monomer ‘m’, solvent ‘s’, or cosolvent ‘c’) and [small sigma, Greek, circumflex]_i is the unit vector representing the orientation of the particle at site i, and is the set of particles with which i interacts. The function ε(·) defines pairwise interaction energies, n is the total number of lattice sites, N_p is the number of polymer chains, N_m is the number of monomers per chain, and r⃑^(k)_i is the position of the lth monomer in chain k. The potential V(·) ensures bonded monomers remain within . The non-bonded interaction term is defined as:where captures the difference between aligned and misaligned interactions, and Λ(i,j) ∈ [0,1] distinguishes the degree of alignment. For convenience, we assume aligned interactions are stronger, so Δ_αiαj ≤ 0. When Λ(i,j) or Δ_αiαj = 0, interactions are isotropic, and FHP reduces to standard FH theory.

We investigate Λ(i,j) with a form given by

Λ_corr(i,j) = Θ([small sigma, Greek, circumflex]_i·[small sigma, Greek, circumflex]_j − δ) (3)

where Θ(·) is the Heaviside step function, and δ defines the angular tolerance for alignment. We refer to this interaction style as a “correlation network,” and it resembles interactions in Maier–Saupe theory.⁷⁵ Correlation networks allow all neighbors to align cooperatively, since the function becomes nonzero when two vectors point nominally in the same direction. Importantly, the orientation vectors here are coarse-grained internal state variables that modulate local interaction energies. In one manifestation, these vectors could encode differences in molecular orientation, where the positioning of certain functional groups would influence the manner of interaction with other moieties. However, these vectors may also represent other forms of collective ordering, such as that seen in nematic liquid crystals or in hydrogen-bond networks in water, where many simultaneous, orientation-dependent interactions occur. While Δ_αiαj terms are not treated in chemically specific terms here, surveying its influence over different values phenomenologically probes the role of orientation-dependent interactions, irrespective of their origin.

We note that the transition from simple FH to FHP introduces several additional parameters governing particle interactions. However, not all parameters are essential for modeling specific phenomena, and FHP parameters have been successfully fitted to experimental data with high precision.⁷⁴ Regularization techniques and physical constraints can yield parsimonious models when needed. Since our objective is to distinguish cononsolvency manifestations, we employ minimal models that vary only two variables while keeping others fixed. In particular, alignment-biased correlation network interactions are assigned between monomer–cosolvent and solvent–cosolvent pairs with fixed weights p_v = 1.0 and p_Ω = 0.25 and the energetic penalty for misalignment is fixed Δ_mc = Δ_sc = Δ; all other interactions are isotropic. This minimal extension is sufficient to capture the anisotropic mechanisms discussed below while avoiding an unwieldy parameter space. To provide some physical context, in the main text, we report results using Δ = 0.5ε₀ where ε₀ = k_BT is the unit of energy. At 300 K, representative non-covalent interactions span a wide energy range: water-to-amide hydrogen bonds are roughly −3 to −10 kcal mol⁻¹,⁷⁶ π–π contacts lie between −3 and −5 kcal mol⁻¹,^77,78 whereas water–alcohol hydrogen bonds are much weaker, about −0.5 kcal mol⁻¹.^79,80 To bracket this spectrum, additional results for Δ = −0.2ε₀, −0.5ε₀, −0.8ε₀ are provided in the SI.

2.2 Helmholtz energy by mean-field analysis

Macroscopic phase behavior in the FHP framework is evaluated in the context of mean-field theory. All systems are considered incompressible, such that the total volume is given by where n_i is the mole number and v_i the molar volume of species i. If each polymer segment occupies one lattice site of volume v₀, then v_p = N_mv₀. The system composition is described by volume fractions ϕ_i = (n_iv_i)/V, with .

Following ref. 74, but extended to three components, the intensive free energy of mixing per lattice site is:

where β = (k_BT)⁻¹ and χ^FHP_ij are FHP interaction parameters (with χ^FHP_ii = 0). This expression mirrors the ternary FH form but with modified interaction terms:where the first term is essentially the FH χ at the misaligned energy scale, z is the coordination number on the lattice, and [small chi, Greek, tilde]_ij(T) is a perturbation term that accounts for the difference between aligned and misaligned interactions. In particular,where [capital Delta, Greek, tilde]_ij(T) accounts for the free energy difference between aligned and misaligned states:

The terms p_v and p_Ω are geometric factors related to what defines aligned interactions. In particular, p_v is the fraction of neighbors that can form aligned interactions, and p_Ω is the fraction of pairwise orientations that are classified as aligned for such neighbors. While eqn (4) retains the classical FH functional form dictated by symmetry and thermodynamic consistency, eqn (5) and (6) allow χ-values to incorporate state-dependent interactions into the FH framework without modifying its structure. Thus, the key offering of leveraging the FHP conceptual framework is that the various terms that contribute to χ^FHP_jk are traceable to a well-defined and immutable microscopic Hamiltonian with controllable orientation-dependent energy terms.

For simplicity, all numerical results correspond to a simple cubic lattice where includes the nearest, next-nearest, and next-next-nearest neighbors of site i, yielding z = 26. Furthermore, polymers are monodisperse with N_m = 72, such that the volume occupied by a polymer chain is v_p = N_mv₀. The volume occupied by each solvent particle and cosolvent particle is v₀. Particle orientations are restricted to the 26 lattice directions. Although choices regarding the lattice and the interaction neighborhood influences z, a different prescription is expected to only modify the energy scales at which phenomena are observed without qualitatively altering mechanisms.

2.3 Characterization of phase behavior

The phase behavior of the system can be determined from eqn (4) by specifying the composition ϕ = (ϕ_p,ϕ_s,ϕ_c), the molar volumes v = (v_p,v_s,v_c), and the interaction parameters χ^FHP = (χ^FHP_ps, χ^FHP_pc, χ^FHP_sc). Briefly, the stability of a homogeneous mixture is assessed by examining the sign of the determinant of the Hessian for the Helmholtz energy:The spinodal boundary is defined by the locus of points for which |H| = 0, while negative values indicate an unstable mixture that would undergo phase separation. Here, we use the condition |H| < 0 over some composition range as the criterion for determining whether a system exhibits cononsolvency at a macroscopic level.

For select systems that meet the criterion above, their phase behavior is more precisely characterized by mapping binodal boundaries and determining the compositions of coexisting phases. Binodal boundaries are computed by first identifying critical points using constraints involving third-order derivatives:^81,82

andFollowing identification of a critical point, coexistence lines are extended using a recently proposed natural parameter continuation algorithm⁸³ akin to Gibbs–Duhem integration. As this approach iteratively identifies different compositions with equivalent chemical potentials, it also facilitates the construction of tie lines to indicate the equilibrium compositions of the coexisting phases. For additional details regarding this algorithm or the utilization of eqn (9)–(11), readers are referred to ref. 83 and the references therein.

2.4 Monte Carlo simulations

To characterize cononsolvency at a microscopic level, we perform lattice Monte Carlo (MC) simulations. The MC simulations directly implement specific manifestations of the Hamiltonian defined by eqn (1). These simulations are performed only for systems with parameters that exhibit mixture instability (Table 1). All simulations feature a 34 × 34 × 34 simple cubic lattice with periodic boundaries and unit edge length; lattice sites are occupied by a mixture of solvent and cosolvent particles and a single polymer chain (N_m = 32). For each manifestation of a Hamiltonian, simulations are performed over a set of cosolvent fractions that approximately (within the limits afforded by the discrete lattice) correspond to x_c ∈ [0,1] in increments of 0.1.

Table 1 Parameters used in Monte Carlo simulations of FHP models. Parameters for label type are for isotropic simulations (Δ_ij = 0) and those for type R_ij are for anisotropic simulations (Δ_ij ≠ 0)

Label	χ^FHPps	χ^FHPpc	χ^FHPsc	Δmc	Δsc
	−1	−2	−10	0	0
R^⋏sc	−1	−2	−10	−0.5	−0.5
	−1	−10	0	0	0
R^⋏pc	−1	−10	0	−0.5	−0.5

---

Systems are initialized by placing polymer chains on lattice sites using a self-avoiding random walk. Cosolvent particles are then randomly distributed across the remaining sites to achieve the target mole fraction, and any unoccupied sites are filled with solvent particles. All particles are randomly assigned one of twenty six possible orientations with uniform probability. Configurational sampling is performed using a variety of MC moves. For solvent and cosolvent particles, these include orientation exchanges, collective orientation perturbations, and particle swaps. Polymer moves include end-rotation, forward and backward reputation, and chain regrowth with orientation updates. Orientation updates are also applied to solvent particles in contact with the polymer and to a randomly selected subset of lattice particles. See Table S1 in the SI for a summary of all moves. All move types are attempted with equal probability.

Each simulation consists of 10⁸ MC moves, with configurations sampled every 10⁴ moves. The first half of each simulation is used for equilibration and the second half for collecting production data (see SI, Fig. S3 and S4). For each condition (parameter set, composition, and temperature), thirty independent simulations are performed for statistical analysis and uncertainty quantification. In addition, to compute certain properties for an ideal mixture, thirty independent simulations are also run under athermal conditions (i.e., all energy parameters except the bonding and implied excluded-volume interactions are set to zero).

It is worth nothing that the MC simulations technically reflect results for a chain in the canonical ensemble. Consequently, it is conceivable that substantial partitioning of solvent species into the near-space volume occupied by the polymer could reduce the effective (co)solvent mole fractions in the “bulk,” and such behavior would depend on the finite size of the simulation cell. Here, we verified that the number of solvent and cosolvent molecules coordinated within the polymer solvation shell remains small relative to the total particle count (<0.5%, see SI, Fig. S2), ensuring that bulk composition remains effectively constant, and this should be sufficient to effectively draw the connection between macroscopic thermodynamics and the dilute single-chain physics. In the future, it may be preferable to employ a constant chemical potential ensemble,^14,45 where the polymer gyration volume exchanges solvent and cosolvent with an external reservoir, though this may require additional considerations to be compatible with the Hamiltonians explored here.

2.5 Conformational analysis

Polymer size is characterized by the radius of gyration R_g. For a polymer with N_m monomers at positions r_i, the R_g iswith the polymer center of mass defined per usualFor consistent comparisons across conditions, we report a normalized radius of gyration, , where 〈R_g〉 is the ensemble-averaged R_g for a given system and 〈R_g〉_ath is a reference given by an athermal simulation incorporating only excluded-volume interactions. For chains with N_m = 72, we empirically find that 〈R_g〉_ath = 6.15 ± 0.005. For simplicity, we use 〈R_g〉_ath = 6.15 in out calculations.

2.6 Solvation analysis

To distinguish cononsolvency mechanisms at the single-chain level, we compare the local solvation environment of the polymer from Monte Carlo simulations with ideal-mixing predictions. Deviations are quantified by the normalized excess number of interactions

N^ex_κ = N_κ − N^id_κ (14)

where N_κ specifies the number of interactions under some restriction κ, and the superscripts ‘ex’ and ‘id’ indicate excess and ideal values, respectively. The restriction of κ for our analysis relates to some pairing of specific species or particles and/or whether interactions are aligned versus misaligned. In the following, numbers of unique pairwise interactions between monomer, solvent, and cosolvent particles are represented with N_αiαj where α_i and α_j are the chemical species p, s, and c.

For ideal references, we first empirically determine N^id_mm from MC simulation of an athermal chain. Subsequently, the quantities N^id_ms and Nⁱ_mc are inferred via

N^id_ms = (1 − x_c)N^id_ms∪mc, (16)

N^id_mc = x_cN^id_ms∪mc, (17)

and

N^id_ms∪mc = zN_m − 2N^id_mm, (18)

where the last expression for the number of either monomer–solvent or monomer–cosolvent interactions exploits the restriction that the overall number of interactions involving the polymer is conserved on the lattice. For chains with N_m = 72, we estimate N^id_mm by computing the average number of monomer–monomer interactions from athermal simulations of a chain in solvent. These simulations yield 〈N_mm〉_ath = 150.12 ± 0.05. For simplicity, we use N^id_mm = 150 in all calculations. Similarly, the expected number of solvent–cosolvent interactions is approximated as

N^id_sc = zN_sx_c. (19)

To assess the structure of the solvation shell, we compute the fraction of ‘bridging’ cosolvents Φ^⋏_b, ‘mediating’ cosolvent Φ^⋏_m, and ‘passive’ osolvents Φ^⋏_p in the solvation shell. The solvation shell is defined as the set of solvent and cosolvent particles in direct contact with the polymer i.e., (nearest, next-nearest, or next-next-nearest neighbor). The classification of any given solvent as bridging, mediating, or passive depends on its underlying Hamiltonian. For systems with orientation-dependent interactions, a bridging cosolvent is aligned and in contact with, two monomer units that lie at least three bonds apart, thereby forming an enthalpic bridge between distant chain segments (Fig. 1(a)). A mediating cosolvent aligns with at least one monomer and one solvent particle but does not meet the bridging criterion. All other cosolvents are labeled passive (Fig. 1(b)). For systems with only isotropic interactions, we make analogous classifications based on the proximity and composition of the solvation shell, but the conditions regarding orientation are relaxed. In other words, a bridging cosolvent contacts two monomers separated by at least three bonds; a mediating cosolvent contacts one monomer and one solvent particle without bridging; and the remainder are passive . After classifying the various cosolvents, we report the normalized mean energy per particle for bridging and mediating cosolvents, denoted Ē_b and Ē_m, respectively; the passive set provides a baseline for comparison.

Fig. 1 Visualization of bridging and mediating cosolvents (orange) interacting with solvents (blue) and a polymer (green) in anisotropic (⋏) systems from Monte Carlo simulations. Monomer segments darken in shade along the contour of the polymer backbone. (a) Bridging cosolvents coordinating and aligning with distal monomer units (shown by gold bonds). (b) Mediating cosolvents coordinating aligned interactions with solvent particles (shown by red bonds) and with monomer units (shown by green-yellow bonds).

3 Results and discussion

Our analysis proceeds in two parts. In Section 3.1, we use mean-field theory to identify the energetic regimes that drive cononsolvency through macroscopic phase separation. These results establish conceptual scaffolding and terminology consistent with prior literature. While similar results have appeared variously,^{34,50–52,68,71} we include this analysis for completeness since application of FHP is new, and specific results here directly inform our subsequent system selection. Readers familiar with cononsolvency may skim this section. In Section 3.2, we employ lattice Monte Carlo simulations to examine microscopic behavior. We strategically select two pairs of systems with equivalent χ^eff_ij from phase-separating regimes, but for a given pair, system differ in their Hamiltonians; one features only isotropic interactions and the other incorporates orientation-dependent terms. Through detailed analysis of coil–globule–coil transitions, which are observed for all systems, we elucidate how interaction anisotropy influences cononsolvency at the microscopic level.

3.1 Analysis of phase behavior

To establish essential conceptual scaffolding and expectations surrounding cononsolvency, we begin by identifying parameter regimes that display phase separation in ternary polymer solutions in the context of the FHP framework. In particular, we systematically explore the interaction-parameter space χ^FHP= (χ^FHP_ps, χ^FHP_pc, χ^FHP_sc) and identify parameter sets for which |H| < 0, signaling instability and thus phase separation, under the constraint that all components are mutually miscible (χ^FHP_ij < 0 ∀i, j).

Consistent with prior literature using FH theory,^{34,50,51,68,71} there are three distinct parameter regimes characterized by the strongest interaction type present (Fig. 2). We denote these regimes as R_ij, corresponding to the region of parameter space where (χ^FHP_ij < χ_ik^FHP, χ_jk^FHP ≤ 0) (i.e., the strongest attractive interaction is between species i and j). The emergence of cononsolvency, in terms of competing interactions, can be understood by examining successive cross sections in Fig. 2(a). Initially, at χ^FHP_sc = 0, two distinct regions, R_ps and R_pc, exhibit phase separation for a broad range of parameters, for which polymer–solvent or polymer–cosolvent interactions are strongly favorable relative to interactions between solvent and cosolvent (Fig. 2(b)). As solvent–cosolvent interactions become increasingly favorable (more negative χ^FHP_sc), the areas of R_ps and R_pc shrink, as interactions are more balanced. Concurrently, a third region, R_sc, emerges where solvent–cosolvent interactions dominate. (Fig. 2(c)) Eventually, as solvent–cosolvent affinity strengthens further, R_ps and R_pc vanish entirely, leaving an expanded R_sc region (Fig. 2(d) and (e)). This analysis provides a baseline expectation that cononsolvency can be driven by different dominant interactions, be it polymer–solvent, polymer–cosolvent, or solvent–cosolvent, irrespective of microscopic details.

Fig. 2 Classification of phase behavior for ternary polymer solutions described by the FHP framework. (a) Regions in interaction parameter space χ^FHP = (χ^FHP_ps, χ^FHP_pc, χ^FHP_sc) where phase separation is required at some composition ϕ, as determined by the mean-field stability criterion (|H| < 0). Markers indicate parameter sets that yield instability; cross-sectional planes at χ^FHP_sc= 0, −4, −7, and −10 correspond to panels (b)–(e). (b)–(e) Two-dimensional slices of (a) showing phase behavior as a function of χ^FHP_ps and χ^FHP_pc for fixed values of χ^FHP_sc. Shaded regions indicate parameter combinations that guarantee phase separation at some composition. Regions of observed two-phase coexistence are labeled R_sc, R_ps, and R_pc, according to the most negative interaction parameter. Square and diamond markers in (b) and (e) denote parameter sets selected for further analysis in Section 3.2.

Prior work has suggested that the shape, location, and extent of coexistence regions, as well as the orientation of tie lines, may reflect distinct underlying mechanisms of cononsolvency.⁵¹ To investigate this, we constructed full ternary phase diagrams for representative parameter sets drawn from each of the three regimes: R_ps, R_pc, and R_sc. Indeed, the shapes of binodal curves and the orientations of tie lines differ based on the dominant interaction present in the system (SI, Fig. S1). When interactions between a polymer and one of the solvents dominate, tie lines connect phases in which one is rich in polymer and the better solvent, while the other is rich in the lesser solvent and depleted in polymer. When solvent–cosolvent interactions dominate, the tie lines orient nearly perpendicular to the polymer composition axis, leading to coexistence between a very polymer-lean phase with mixed solvent and a polymer-rich phase.

Visually, the orientation of the tielines and the resultant coexistent composition suggest that there is a polymer-lean and polymer-rich phase. To verify that the reported instabilities correspond to cononsolvency, we performed an eigenmode analysis of the Hessian of the free energy at representative compositions within the binodal region. In the χ_pc-dominated regime, the unstable mode corresponds to a channel where the polymer and cosolvent fluctuate together while the solvent fluctuates in the opposite direction. Conversely, in the χ_sc-dominated regime, the instability primarily follows the polymer concentration while solvent and cosolvent fluctuate in phase. These results confirm that cononsolvency arises via distinct mechanisms in the two regimes. Full details of the eigenmode decomposition and channel classification are provided in the SI (Section S2).

Nevertheless, while these observations may serve as useful indicators of the dominant underlying interactions, they remain qualitative and do not offer insight into the microscopic details of the interactions that underlie cononsolvency.

3.2 Analysis of simulations

The preceding results demonstrate how specific combinations of mean-field interaction parameters can induce macroscopic phase separation across a range of ternary compositions. We now turn to the question of whether such macroscopic behavior is echoed by a coil–globule transition induced by solvent-mixture composition and how this depends on whether interactions are isotropic or anisotropic in form.

To address this, we perform Monte Carlo (MC) simulations to study the conformational behavior and solvation-shell characteristics of a single polymer chain in mixed solvents for four systems. We present results for a chain with degree of polymerization N_m = 72; corresponding data for a shorter chain with N_m = 32 are provided in the SI (Fig. S9–S12). Both lengths show similar physics but the longer chain shows more pronounced effects. In Section 3.2.1, we describe the selection of these systems and broadly categorize their behavior. The results supporting this categorization are then detailed in Sections 3.2.2 and 3.2.3.

3.2.1 Overview of systems and mechanisms. The four systems studied divide into two pairs based on their parameter sets (Table 1). One pair shares the same χ^FHP_ij parameter set from R_sc (diamond marker in Fig. 2(e)), and the other shares χ^FHP_ij from R_pc (square marker in Fig. 2(b)). For a given pair, systems differ in whether non-bonded interactions are orientation-independent (isotropic, denoted with ‘°’) or orientation-dependent (anisotropic, denoted with ‘⋏’). The complexity of anisotropy is restricted to only interactions involving cosolvent. For simplicity, we impose that this is governed by a single energy scale, Δ_*c = Δ_mc = Δ_sc. In the main text, we focus on Δ_*c = 0.5, while results for lesser and greater Δ_*c are provided in the SI. Furthermore, parameters associated with monomer–monomer, solvent–solvent, cosolvent–cosolvent and monomer–solvent interactions are fixed across all systems (, , and ).

In all systems, we find evidence of single-chain collapse at certain mixture compositions based on observed R_g relative to that in either pure solvent. Based on further analysis of the solvation environment around the polymer, we then associate these chain collapses with specific mechanisms within the vernacular of the cononsolvency literature. Simulations with parameter sets that are derived from R_sc show a collapse through preferential mixing of solvent and cosolvent; this mechanism tends to yield a “dry” globule. In this case, including anisotropic interactions with cosolvent species qualitatively alters the cosolvent-induced coil–globule transition. Parameter sets from R_pc collapse through preferential adsorption of cosolvent; this mechanism tends to yield a “wet” globule. For this situation, the role of anisotropy is apparent only in nuanced changes to the polymer solvation environment in its collapsed state.

3.2.2 Chain collapse by preferential mixing of solvent and cosolvent. Chain collapse from an extended solvated polymer to a dry, solvent-excluded globule is observed in systems where solvent–cosolvent interactions are strongly favorable (R_sc). In the isotropic case (, blue circles), the collapse coincides with nearly equimolar mixtures at which entropy of mixing is maximized, resulting in a symmetric R_g profile as a function of cosolvent mole fraction (Fig. 3(a)); the collapse is likewise evident by the enrichment in monomer–monomer interactions (Fig. 3(b)), which is also nearly symmetric. The resulting globule is dry, as indicated by the depletion in excess monomer–cosolvent interactions (Fig. 3(c)). The excess number of monomer–cosolvent contacts remains small, indicating behavior close to that expected from ideal mixing and thus limited preferential interaction between cosolvent and polymer. Furthermore, the chain is maximally collapsed at x_c ≈ 0.5, which aligns with the excess number of solvent–cosolvent interactions (Fig. 3(d)). This dry globule and congruence of polymer collapse with enrichment in solvent–cosolvent interactions are defining characteristics of cononsolvency driven by solvent–cosolvent interactions.

Fig. 3 A comparison of conformational characteristics and environment of a polymer from (blue circles) and R_sc (orange triangles). (a) Normalized single-chain R_g and normalized excess (b) monomer–monomer, (c) monomer–cosolvent and (d) solvent–cosolvent interactions as cosolvent fraction x_c is varied from 0 to 1. In (a), the dashed horizontal line is a guide for the reference to a maximally compact polymer and the white arrow shows the shift in the minima. In (a) and (b), the white arrow shows a marked shift in the behavior of the curves. Horizontal axis labels are share between panels (a) and (c) as well as (b) and (d). Error bars correspond to the standard error of the mean and are generally smaller than the symbol size.

However, there are notable differences in the anisotropic case (R_sc, orange triangles). In particular, the maximal chain collapse shifts to lesser cosolvent fractions—in this case, around x_c ≈ 0.3 (Fig. 3(a)). This is consistent with the maximum in the number of monomer–monomer contacts (Fig. 3(b)), but it is inconsistent with maximum enrichment of solvent–cosolvent interactions, which remains at equimolar concentrations. The collapsed state of the polymer remains dry, largely excluding both solvent and cosolvent. Collectively, these observations highlight that the distinguishing feature between these two systems is the composition at which the polymer is maximally collapsed. Whereas both systems should effectively share the same macroscopic phase behavior (at the given temperature), systems with strong anisotropy will be characterized by a more asymmetric R_g profile as a function of cosolvent mole fraction.

To better understand what drives this asymmetric collapse, we perform a configurational and energetic analysis of the polymer solvation shell. In both systems, on trace addition of cosolvent, nearly all cosolvent particles in the solvation shell act as mediating cosolvents, coordinating interactions between monomer and solvent particles, with the remainder acting as bridging cosolvents (Fig. 4(a) and (b)). The relative populations of cosolvent particles and how they interact with the polymer is thus not particularly distinctive. However, the relative energetic contributions from these groups does distinguish the two systems effectively. Each bridging and mediating cosolvent in the anisotropic system provides substantially stronger stabilizing interactions to the polymer over certain composition ranges (Fig. 4(c) and (d)). The difference in stabilization between isotropic and anisotropic systems is largest around x_c ∈ [0.2,0.4], thereby providing a substantial enthalpic driving force at lower cosolvent fractions and shifting the collapse from x_c = 0.5 in the isotropic case to x_c = 0.3 in the anisotropic case.

Fig. 4 Analysis and comparison of the solvation environment of systems and R^⋏_sc. (a) Fraction of bridging cosolvents Φ_b, (b) fraction of mediating cosolvents Φ_m, (c) energetic contribution per cosolvent particle in the solvation shell Ē_b, and (d) the energetic contribution per cosolvent particle in the solvation shell Ē_m as cosolvent fraction x_c is varied from 0 to 1. Horizontal axis labels are share between panels (a) and (c) as well as (b) and (d). The white arrows highlight a marked difference in the magnitude of intensive energies. Error bars correspond to the standard error of the mean and are generally smaller than the symbol size.

Ultimately, this elucidates the qualitative difference between single-chain collapse in systems where solvent–cosolvent interactions dominate. The correlation network formed by bridging and mediating interactions stabilizes the polymer sufficiently to overcome configurational entropy loss, yielding asymmetric collapse. By contrast, isotropic systems require maximized solvent mixing to compensate for reduced polymer configurational entropy.

3.2.3 Chain collapse by preferential adsorption of cosolvent. In a regime dominated by polymer–cosolvent affinity (R_pc), most signatures of microscopic cononsolvency are effectively equivalent between systems with equivalent χ^FHP_ij, irrespective of whether such parameters arise from purely isotropic or include anisotropic terms. For example, both isotropic (, blue circles) and anisotropic (R^⋏_pc orange triangles) systems display asymmetric cosolvent-induced coil–globule transitions, as monitored by R_g (Fig. 5(a)). Notably, this asymmetric profile may be reminiscent of that found for the anisotropic system in Fig. 3(a). Nevertheless, by comparing the excess number of monomer–monomer interactions in Fig. 3(b) versus those in Fig. 5(b), one can infer that the collapsed chain induced by polymer–cosolvent interactions is relatively wet. In particular, the globules are found to be cosolvent-laden globule, as evidenced by an excess of cosolvent–monomer interactions (Fig. 5(c)). Between isotropic and anisotropic systems, this enrichment is slightly stronger in the latter. The simulations with anisotropic interactions display more significant solvent–cosolvent interactions at intermediate cosolvent fractions (Fig. 5(d)). However, this enhanced mixing is of little consequence to the polymer conformational behavior and simply arises to due Δ_sc being negative. Because of the strong enrichment in polymer–cosolvent interactions, we mechanistically refer to behavior in this regime as preferential adsorption.

Fig. 5 A comparison of conformational characteristics and environment of a polymer from (blue circles) and R^⋏_pc (orange triangles). (a) Normalized single-chain R_g and normalized excess (b) monomer–monomer, (c) solvent–cosolvent, and (d) monomer–cosolvent interactions as cosolvent fraction x_c is varied from 0 to 1. In (a), the dashed horizontal line is a guide for the reference to a maximally compact polymer. Horizontal axis labels are share between panels (a) and (c) as well as (b) and (d). In (c) and (d), the white arrow shows a marked shift in the trendlines. Error bars correspond to the standard error of the mean and are generally smaller than the symbol size.

The role of anisotropy for preferential-adsorption driven cononsolvency appears limited to slight differences in solvation motifs. On addition of cosolvent, both systems with isotropic and anisotropic interactions possess significant fractions of bridging and mediating cosolvents (Fig. 6(a) and (b)). By contrast to behavior in the preferential mixing regime, however, here there are notable differences in relative proportions. For the system with purely isotropic interactions, at low cosolvent fractions, the number of mediating cosolvents exceeds the number of bridging solvents. However, making interactions with cosolvent orientation-dependent reverses this trend, such that more cosolvent particles are interacting simultaneously with distal monomers on the polymer chain, rather than being situated between a monomer and solvent. In conjunction with Fig. 5(c), this implies that the globule in the anisotropic system has more cosolvent embedded or intercalated within its pervaded volume. Finally, the energetic contribution per cosolvent particle are also notably enhanced when including anisotropic interactions (Fig. 6(c) and (d)). While this does not seemingly have any clear impact on the presence or nature of the coil–globule transition, it may have implications on the properties of the collapsed state.

Fig. 6 Analysis and comparison of the solvation environment of systems and R^⋏_pc. (a) Fraction of bridging cosolvents Φ_b, (b) fraction of mediating cosolvents Φ_m, (c) energetic contribution per cosolvent particle in the solvation shell Ē_b and (d) the energetic contribution per cosolvent particle in the solvation shell Ē_m as cosolvent fraction x_c is varied from 0 to 1. Horizontal axis labels are share between panels (a) and (c) as well as (b) and (d). In (c) and (d), the white arrow shows a marked shift in the intensive energy. Error bars correspond to the standard error of the mean and are generally smaller than the symbol size.

Thus, cononsolvency via preferential adsorption exhibits qualitatively similar behavior for both isotropic and anisotropic interactions. In both cases, cosolvent particles percolate the gyration volume, expel solvent, and intercalate between monomer segments to form wet globules. This supports prior observations⁵⁶ that a polymer chain can collapse without reduction in overall solvent quality. In this case, this arises because the cosolvent intercalates with strong affinity to the polymer chain. When interactions are only isotropic, this produces abrupt collapse where cosolvent mediates monomer–solvent interactions; the addition of more cosolvent leads to fewer such interactions, causing gradual re-expansion. Anisotropic interactions intensify these physics through enhanced cosolvent intercalation and stronger mediating interactions that persist throughout the transition. While the underlying physical drivers remain effectively identical, structural and compositional correlations in globules induced by polymer–cosolvent affinity may differ for sufficiently strong anisotropic interactions.

4 Conclusions

We investigated how orientational interactions influence cononsolvency phenomena in ternary polymer solutions. Using the Flory–Huggins–Potts (FHP) framework enabled systematic comparison of systems with identical effective χ parameters but different underlying interaction types, including those with orientation-dependent energetic contributions. After identifying phase-separation regimes driven by either solvent–cosolvent or polymer–cosolvent affinity, we contrasted systems with the same effective χ parameters but achieved either purely through isotropic interactions or with inclusion of anisotropic interactions, specifically with cosolvent species. This enabled controlled study of how anisotropic interactions affect microscopic physics despite equivalent macroscopic phase behavior. Analysis focused on cosolvent-induced coil–globule transitions. This work complements prior extensive literature on cononsolvency based on FH theory.

Our analysis revealed that importance of orientation-dependent interactions depends on the dominant interaction type driving cononsolvency. For systems driven by strong solvent–cosolvent interactions, orientation-dependent interactions qualitatively altered coil–globule transitions, inducing collapse at lower cosolvent fractions than isotropic systems with equivalent mean-field parameters. This enhanced, asymmetric collapse was attributed to stronger cosolvent-mediated interactions, where cosolvent acts as a surfactant favorably interacting with both polymer and solvent.⁴⁷ This asymmetric collapse notably resembled that observed in polymer–cosolvent affinity-driven systems. However, whereas solvent–cosolvent affinity produced dry, solvent-depleted globules, those produced by polymer–cosolvent affinity were more wet with embedded cosolvent; this distinction offers a pathway to distinguish these regimes at the microscopic level. When comparing isotropic and anisotropic systems both dominated by polymer–cosolvent affinity, anisotropic interactions primarily increase bridging-type cosolvent configurations and potentially enhance globule stability, without fundamentally altering the mechanism.

This work highlights several implications and opportunities for future inquiry. Discriminating among cononsolvency mechanisms proves challenging when examining only macroscopic phase behavior or microscopic single-chain conformations in isolation. However, complementary analysis at both scales suggests that preferential mixing, anisotropically-influenced preferential mixing, and preferential adsorption mechanisms can be distinguished. While investigating the role of anisotropy interactions via experiment remains nontrivial, our results suggest that temperature-dependent characterization across compositions may be informative. While the FHP framework provides physically grounded parameters enabling precise experimental fitting,⁷⁴ this has so far relied on macroscopic phase behavior. Connecting FHP parameters to molecular simulations or microscopic experimental observables would clarify when orientation-dependent interactions emerge in real systems. Meanwhile, we propose that studies employing systematic molecular modifications—such as adding electron-withdrawing substituents or stackable molecular moieties to cosolvents to tune hydrogen-bonding and molecular packing—may offer tangible routes to probe anisotropic effects without requiring FHP abstraction.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data and code associated with this study are publicly accessible at https://github.com/webbtheosim/MCLATTE-public.

Supplementary information (SI): Monte Carlo simulation details; representative phase diagrams; eigenmode analysis of the Hessian; robustness and convergence checks; details of microscopic models; additional results for other parameter combinations. See DOI: https://doi.org/10.1039/d5sm00767d.

Acknowledgements

S. D. and M. A. W. acknowledge support from the National Science Foundation under grant no. 2237470. Calculations were performed using resources from Princeton Research Computing at Princeton University, which is a consortium led by the Princeton Institute for Computational Science and Engineering (PICSciE) and Office of Information Technology's Research Computing.

Notes and references

1. M. Heskins and J. E. Guillet, J. Macromol. Sci., Part A: Pure Appl. Chem., 1968, 2, 1441–1455.
2. S. Bharadwaj, B.-J. Niebuur, K. Nothdurft, W. Richtering, N. F. A. van der Vegt and C. M. Papadakis, Soft Matter, 2022, 18, 2884–2909.
3. N. Takahashi, T. Kanaya, K. Nishida and K. Kaji, Polymer, 2003, 44, 4075–4078.
4. D. Jia, M. Muthukumar, H. Cheng, C. C. Han and B. Hammouda, Macromolecules, 2017, 50, 7291–7298.
5. I. Fernández-Piérola and A. Horta, Polym. Bull., 1980, 3, 273–278.
6. F. Wang, Y. Shi, S. Luo, Y. Chen and J. Zhao, Macromolecules, 2012, 45, 9196–9204.
7. G. Zhang and C. Wu, J. Am. Chem. Soc., 2001, 123, 1376–1380.
8. Y. Zhao and K. Kremer, Macromol. Rapid Commun., 2022, 43, 2100907.
9. Y. Zhao and K. Kremer, J. Phys. Chem. B, 2021, 125, 9751–9756.
10. J. Hao, H. Cheng, P. Butler, L. Zhang and C. C. Han, J. Chem. Phys., 2010, 132, 154902.
11. S. Sun and P. Wu, Macromolecules, 2010, 43, 9501–9510.
12. M. J. A. Hore, B. Hammouda, Y. Li and H. Cheng, Macromolecules, 2013, 46, 7894–7901.
13. D. Jia, T. Zuo, S. Rogers, H. Cheng, B. Hammouda and C. C. Han, Macromolecules, 2016, 49, 5152–5159.
14. Y. A. Budkov and M. G. Kiselev, J. Phys.: Condens. Matter, 2017, 30, 043001.
15. F. Pooch, V. Teltevskij, E. Karjalainen, H. Tenhu and F. M. Winnik, Macromolecules, 2019, 52, 6361–6368.
16. T. Zuo, C. Ma, G. Jiao, Z. Han, S. Xiao, H. Liang, L. Hong, D. Bowron, A. Soper, C. C. Han and H. Cheng, Macromolecules, 2019, 52, 457–464.
17. H. Yong, E. Bittrich, P. Uhlmann, A. Fery and J.-U. Sommer, Macromolecules, 2019, 52, 6285–6293.
18. Y. Zhang, W. S. Carvalho, C. Fang and M. J. Serpe, Sens. Actuators, B, 2019, 290, 520–526.
19. S. Nian and L. Pu, J. Org. Chem., 2018, 84, 909–913.
20. D. Kleinschmidt, K. Nothdurft, M. V. Anakhov, A. A. Meyer, M. Mork, R. A. Gumerov, I. I. Potemkin, W. Richtering and A. Pich, Mater. Adv., 2020, 1, 2983–2993.
21. Y. Yu, M. Cirelli, B. D. Kieviet, E. S. Kooij, G. J. Vancso and S. de Beer, Polymer, 2016, 102, 372–378.
22. Q. Chen, E. S. Kooij, X. Sui, C. J. Padberg, M. A. Hempenius, P. M. Schön and G. J. Vancso, Soft Matter, 2014, 10, 3134.
23. X. Wang, H. Huang, H. Liu, F. Rehfeldt, X. Wang and K. Zhang, Macromol. Chem. Phys., 2019, 220, 1800562.
24. H. Yong, B. Molcrette, M. Sperling, F. Montel and J.-U. Sommer, Macromolecules, 2021, 54, 4432–4442.
25. C. E. Mills, E. Ding and B. Olsen, Ind. Eng. Chem. Res., 2019, 58, 11698–11709.
26. J. J. Valente, K. S. Verma, M. C. Manning, W. William Wilson and C. S. Henry, Biophys. J., 2005, 89, 4211–4218.
27. T. Xiao, K. H. Gardner and S. R. Sprang, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 11151–11156.
28. G. Reddy, A. N. Muttathukattil and B. Mondal, Curr. Opin. Struct. Biol., 2020, 60, 101–109.
29. D. R. Canchi and A. E. García, Annu. Rev. Phys. Chem., 2013, 64, 273–293.
30. R. Gazi, S. Maity and M. Jana, ACS Omega, 2023, 8, 2832–2843.
31. M. A. C. Stuart, W. T. S. Huck, J. Genzer, M. Müller, C. Ober, M. Stamm, G. B. Sukhorukov, I. Szleifer, V. V. Tsukruk, M. Urban, F. Winnik, S. Zauscher, I. Luzinov and S. Minko, Nat. Mater., 2010, 9, 101–113.
32. D. L. Nettles, A. Chilkoti and L. A. Setton, Adv. Drug Delivery Rev., 2010, 62, 1479–1485.
33. J. C. Rodrí­guez-Cabello, F. J. Arias, M. A. Rodrigo and A. Girotti, Adv. Drug Delivery Rev., 2016, 97, 85–100.
34. J. Dudowicz, K. F. Freed and J. F. Douglas, J. Chem. Phys., 2015, 142, 214906.
35. J. Wang, N. Wang, B. Liu, J. Bai, P. Gong, G. Ru and J. Feng, Phys. Chem. Chem. Phys., 2017, 19, 30097–30106.
36. S. Bharadwaj and N. F. A. van der Vegt, Macromolecules, 2019, 52, 4131–4138.
37. H. G. Schild, M. Muthukumar and D. A. Tirrell, Macromolecules, 1991, 24, 948–952.
38. D. Mukherji, C. M. Marques, T. Stuehn and K. Kremer, J. Chem. Phys., 2015, 142, 114903.
39. F. Tanaka, J. Chem. Phys., 1985, 82, 2466–2471.
40. T. Tanaka, I. Nishio, S.-T. Sun and S. Ueno-Nishio, Science, 1982, 218, 467–469.
41. F. Tanaka, T. Koga and F. M. Winnik, Phys. Rev. Lett., 2008, 101, 028302.
42. F. Tanaka, T. Koga, H. Kojima, N. Xue and F. M. Winnik, Macromolecules, 2011, 44, 2978–2989.
43. C. Dalgicdir and N. F. A. van der Vegt, J. Phys. Chem. B, 2019, 123, 3875–3883.
44. H. A. Pérez-Ramí­rez, C. Haro-Pérez, E. Vázquez-Contreras, J. Klapp, G. Bautista-Carbajal and G. Odriozola, Phys. Chem. Chem. Phys., 2019, 21, 5106–5116.
45. Y. A. Budkov and A. L. Kolesnikov, Polym. Sci., Ser. C, 2018, 60, 148–159.
46. Y. A. Budkov and A. L. Kolesnikov, Soft Matter, 2017, 13, 8362–8367.
47. S. Bharadwaj, D. Nayar, C. Dalgicdir and N. F. A. van der Vegt, Commun. Chem., 2020, 3, 1–7.
48. A. Pica and G. Graziano, Phys. Chem. Chem. Phys., 2016, 18, 25601–25608.
49. C. Dalgicdir, F. Rodrí­guez-Ropero and N. F. A. van der Vegt, J. Phys. Chem. B, 2017, 121, 7741–7748.
50. J. Dudowicz, K. F. Freed and J. F. Douglas, J. Chem. Phys., 2015, 143, 131101.
51. P. Zhang, Macromolecules, 2024, 57, 4298–4311.
52. X. Zhang, J. Zong and D. Meng, Soft Matter, 2020, 16, 7789–7796.
53. W. W. Oliver, W. M. Jacobs and M. A. Webb, J. Phys. Chem. B, 2025, 129, 9551–9565.
54. L. Liu, C. Duan and R. Wang, Macromolecules, 2024, 57, 10694–10703.
55. S. Bharadwaj, D. Nayar, C. Dalgicdir and N. F. A. van der Vegt, J. Chem. Phys., 2021, 154, 134903.
56. D. Mukherji, C. M. Marques and K. Kremer, Nat. Commun., 2014, 5, 1–6.
57. J.-U. Sommer, Macromolecules, 2017, 50, 2219–2228.
58. J.-U. Sommer, Macromolecules, 2018, 51, 3066–3074.
59. D. Mukherji and K. Kremer, Macromolecules, 2013, 46, 9158–9163.
60. D. Mukherji, M. Wagner, M. D. Watson, S. Winzen, T. E. de Oliveira, C. M. Marques and K. Kremer, Soft Matter, 2016, 12, 7995–8003.
61. A. Pica and G. Graziano, Soft Matter, 2017, 13, 7698–7700.
62. N. F. A. van der Vegt and F. Rodrí­guez-Ropero, Soft Matter, 2017, 13, 2289–2291.
63. D. Mukherji, M. Wagner, M. D. Watson, S. Winzen, T. E. de Oliveira, C. M. Marques and K. Kremer, Soft Matter, 2017, 13, 7701–7703.
64. P. de Gennes, Phys. Lett. A, 1972, 38, 339–340.
65. P. D. Gennes, J. Phys., Lett., 1975, 36, 55–57.
66. P. de Gennes, J. Phys., Lett., 1978, 39, 299–301.
67. P. G. de Gennes, Scaling Concepts in Polymer Physics, Cornell University Press, 1979.
68. P. Zhang, Macromolecules, 2025, 58, 2472–2483.
69. P. Zhang, Z. Wang and Z.-G. Wang, Macromolecules, 2022, 56, 153–165.
70. D. Marcato, A. Giacometti, A. Maritan and A. Rosa, Phys. Rev. Mater., 2024, 8, 125601.
71. X. Li, Z. Wang, Z. Wang, Y. Yin, R. Jiang, P. Zhang and B. Li, Soft Matter, 2025, 21, 4858–4868.
72. C. Scherzinger, A. Schwarz, A. Bardow, K. Leonhard and W. Richtering, Curr. Opin. Colloid Interface Sci., 2014, 19, 84–94.
73. H. Yong and J.-U. Sommer, Macromolecules, 2022, 55, 11034–11050.
74. S. Dhamankar and M. A. Webb, ACS Macro Lett., 2024, 13, 818–825.
75. P. van der Schoot, Molecular Theory of Nematic (and Other) Liquid Crystals: An Introduction, Springer International Publishing, 2022.
76. D. A. Dixon, K. D. Dobbs and J. J. Valentini, J. Phys. Chem., 1994, 98, 13435–13439.
77. R. G. Huber, M. A. Margreiter, J. E. Fuchs, S. von Grafenstein, C. S. Tautermann, K. R. Liedl and T. Fox, J. Chem. Inf. Model., 2014, 54, 1371–1379.
78. K. Molčanov and B. Kojić-Prodić, IUCrJ, 2019, 6, 156–166.
79. I. A. Finneran, P. B. Carroll, M. A. Allodi and G. A. Blake, Phys. Chem. Chem. Phys., 2015, 17, 24210–24214.
80. P. V. Swathi, U. Abdulkareem, T. R. Kartha and V. Madhurima, ChemistrySelect, 2022, 7, e202200413.
81. R. L. Scott, J. Chem. Phys., 1949, 17, 268–279.
82. H. Tompa, J. Chem. Phys., 1949, 17, 1006.
83. S. Dhamankar, S. Jiang and M. A. Webb, Mol. Syst. Des. Eng., 2025, 10, 89–101.

---