PhaseTime: a Flory–Huggins-based computational framework for predicting time-dependent phase diagrams of reactive polymer blends

Abstract

Reactive polymer blends often undergo phase separation while their molecular constitution is still evolving. In such systems, stability limits and coexistence boundaries are not fixed properties of the initial mixture, but change with reaction progress, temperature, conversion, or another imposed protocol variable. This makes the interpretation and prediction of reaction-induced phase separation challenging, particularly when the components are polydisperse and the coexisting phases may differ in both composition and molecular distribution. Here we present PhaseTime, a Python framework for computing time- or protocol-dependent phase diagrams of reactive polydisperse polymer blends within Flory–Huggins-type thermodynamics. The framework combines evolving molecular-distribution models with calculations of spinodal curves, critical points, binodal curves, and cloud/shadow curves. It supports temperature- and composition-dependent interaction parameters, monodisperse reaction-dependent approximations, and polydisperse Flory–Stockmayer distributions. The framework also includes parameter fitting, diagnostic, and plotting tools through both a command-line interface and a Python API. Representative calculations demonstrate idealized phase-diagram topologies, including UCST, LCST, combined UCST/LCST, hourglass, and closed-loop behavior, as well as fitting to literature data and reaction-dependent phase behavior in monodisperse and polydisperse systems. The resulting phase diagrams provide quasi-equilibrium thermodynamic references for interpreting reaction-induced phase separation and for comparison with spatially resolved models in the fast-demixing limit.

---

Design, System, Application

PhaseTime is a Python-based computational framework that supports the design and interpretation of polymer blend formulations and processing protocols in which phase separation controls final morphology and properties such as toughness, permeability, adhesion, optical/electrical response, heat transport and flow. Instead of computing a single static phase diagram, PhaseTime follows how phase behavior evolves along a prescribed coordinate, such as time, temperature, conversion, or another processing variable. This enables formulation variables, reaction pathways, interaction models, and molecular weight distributions to be examined systematically as design parameters. The current implementation considers quasibinary blends of two homopolymeric species within a Flory–Huggins-type mean-field framework. Each species may be monodisperse or represented by a conversion-dependent molecular-size distribution, including pregelation Flory–Stockmayer distributions. At each protocol state, PhaseTime updates the molecular constitution and thermodynamic parameters, then evaluates spinodal curves, critical points, binodal curves, and cloud/shadow curves. PhaseTime generates thermodynamic design maps for reaction-induced or protocol-driven phase separation. These maps identify when formulations enter metastable or unstable regions, how coexistence compositions evolve during reaction or cooling, and how polydispersity or molecular growth modifies the accessible processing window. The resulting phase diagrams are relevant to modified thermosets, toughened polymer networks, coatings, adhesives, and membranes, and provide reference data for experiments and spatially resolved simulations.

1. Introduction

Polymer blend phase behavior is a central topic in polymer physics, physical chemistry, and materials science. Polymer blending is widely used as a comparatively simple and cost-effective route for modifying material properties for specific applications.^1–4 The miscibility of the components determines whether a mixture remains homogeneous or separates into distinct phases, and therefore influences morphology, interfacial structure, and macroscopic properties.⁴ In thermosetting systems, the incorporation of dispersed modifiers, such as thermoplastics or rubbery components, can improve properties including toughness and fracture resistance.^2,3 Since the morphology is often generated during processing, predicting when and how phase separation occurs is important for relating formulation and processing conditions to final material properties. Phase diagrams provide a useful framework for identifying stability limits, coexistence regions, and critical conditions in multicomponent polymer systems.^4–8 Such information is relevant in the analysis of polymer blends,^4–7,9 rubber-modified thermosets,^2–4,10 membranes,¹¹ hyperbranched polymers,^8,12 adhesives,^13,14 coatings,¹⁵ and more.

Knowledge of stability limits is particularly important in systems undergoing reaction-induced phase separation.^{2–4,8,10,16} In modified thermosetting polymers, for example, an initially homogeneous mixture may become unstable because polymerization or curing changes the molecular-weight distributions and effective interaction parameters of the system.^3,4 Once the reaction path crosses the coexistence boundary, phase separation competes with the continuing reaction. The resulting morphology is therefore controlled by both the evolving thermodynamic phase diagram and the rate of demixing relative to the reaction kinetics.

If demixing is fast compared with reaction, the system may remain close to the quasi-equilibrium coexistence state at each reaction state, so that the phase diagram provides the corresponding coexisting compositions.⁴ Slower demixing allows the system to penetrate into the metastable region before phase separation begins, favoring nucleation and growth.^4,17,18 If the reaction proceeds beyond the spinodal before appreciable demixing occurs, composition fluctuations may grow spontaneously through spinodal decomposition;^4,17–19 the resulting domains may subsequently coarsen through coalescence and Ostwald ripening.^20,21

Near-critical mixtures are more likely to develop co-continuous structures,²² whereas off-critical mixtures more commonly form droplet–matrix morphologies in which the lower-volume phase is usually dispersed within the higher-volume phase. Phase inversion²³ may occur when changes in the coexisting-phase volume fractions cause the phase associated with the initially less abundant component to become continuous. Gelation or vitrification can further arrest demixing and preserve such non-equilibrium morphologies on the processing time scale.^22,24,25 For example, Li et al.²² demonstrated that vitrification of urea-rich domains preserves a co-continuous morphology in flexible polyurethane foams.

In reactive blends, these morphological outcomes are connected to the fact that the thermodynamic phase diagram itself changes during the process. As reaction proceeds, the molecular constitution, molecular-weight distribution, interaction parameters, and phase-equilibrium conditions may evolve, so that the system follows a path through a sequence of instantaneous thermodynamic states. It is therefore useful to consider phase diagrams as functions of a reaction or processing coordinate, such as conversion, temperature, time, or another imposed protocol variable.

From this perspective, reaction- or protocol-dependent phase diagrams provide a thermodynamic reference for interpreting kinetic models of morphology formation. They identify the points along a reaction or processing path where the homogeneous state becomes metastable or unstable, and they provide the corresponding equilibrium or quasi-equilibrium coexistence information. In polydisperse systems, this information may also include changes in the molecular distributions of the coexisting phases, since fractionation can affect coexistence curves and critical behavior.

Flory–Huggins-type thermodynamics^26,27 provides a common theoretical basis for describing polymer-blend phase behavior. In its classical form, the model combines a combinatorial entropy of mixing with an effective interaction parameter and offers a simple mean-field route for calculating stability limits, critical points, and coexistence curves.⁴ These calculations are relatively direct for monodisperse binary blends, where the thermodynamic state is described by a single composition variable. Analytical solutions have been proposed in recent years for two-component²⁸ and multicomponent systems.²⁹ They become more complex for polydisperse systems, where each polymeric species is represented by a distribution of molecular constituents. In such systems, phase separation can involve fractionation between coexisting phases, leading to distinct cloud and shadow curves rather than a single binary binodal. Cloud- and shadow-curve calculations for polydisperse polymer systems have been developed through several classical approaches, including the method of Kamide and co-workers,^30–32 Koningsveld and Staverman,³³ Šolc,³⁴ and Mumby et al.⁶ These methods share the same thermodynamic basis: the cloud point is the state of a prescribed parent phase at which an infinitesimal shadow phase can coexist with it.

Existing open-access computational tools cover related parts of the problem. Web-based implementations, such as the 3PDB Flory–Huggins application³⁵ and interactive demonstrations,³⁶ provide accessible calculations of binary Flory–Huggins free-energy curves and phase diagrams, including spinodal and binodal information. More general open-source packages have also appeared; for example, Flory³⁷ provides numerical tools for finding coexisting phases in multicomponent mixtures with Flory–Huggins-type free energies and related generalizations. These tools are useful for standard or generalized Flory–Huggins calculations, but they are not primarily designed for reactive or protocol-dependent polydisperse blends, and they do not provide the classical cloud- and shadow-curve calculations used for quasibinary polydisperse polymer systems. In such systems, molecular distributions, interaction parameters, stability limits, and cloud/shadow curves must be updated consistently along a reaction or processing path. This creates a more specific computational workflow, in which phase-equilibrium calculations must be coupled to evolving molecular constitution and repeated over a prescribed protocol coordinate.

The present work develops PhaseTime to address this workflow. The framework combines evolving molecular-distribution models with Flory–Huggins-type phase-equilibrium calculations, within a single Python implementation, with the cloud/shadow construction following the formulation of Mumby and coworkers.^5,6 The protocol coordinate may represent time, temperature, conversion, or another externally prescribed variable, allowing stability and coexistence information to be followed along a reaction or processing path. User-defined expressions may be supplied for the evolution of temperature, T(t), and for the extents of reaction of the two species, a₁(t) and a₂(t). At each protocol state, PhaseTime updates the molecular constitution and thermodynamic parameters, and computes spinodal, critical, binodal, cloud, and shadow curves.

The framework is modular with respect to the main model components. The interaction parameter may include separate temperature- and composition-dependent contributions, and additional functional forms can be introduced within the same structure. The implemented temperature-dependent form³⁸ can reproduce several common phase-diagram topologies, including UCST, LCST, combined UCST/LCST behavior, closed-loop, and hourglass-type diagrams. For strongly composition-dependent interactions,³⁸ the framework can also identify multiple critical points and hidden coexistence branches. Molecular constitution may be represented by monodisperse reaction-dependent expressions, such as a generalized Carothers-type relation,^3,4,39 or by polydisperse distributions, including the general Flory–Stockmayer distribution.^40,41

The computed phase diagrams represent the quasi-equilibrium limiting case in which demixing relaxes much faster than the imposed protocol. PhaseTime does not address spatial morphology evolution directly, but describes the local thermodynamic equilibrium problem underlying phase-field or Cahn–Hilliard-type models in the rapid-relaxation limit.^42–44 This makes the results useful as reference phase diagrams for interpreting experiments, benchmarking thermodynamic assumptions, and comparing with spatially resolved simulations of reaction-induced phase separation.^42–46

The remainder of this article is organized as follows. Section 2 presents the theoretical and numerical framework, including the quasibinary thermodynamic formulation, Flory–Huggins free energy, interaction models, molecular-distribution functions, phase-stability and coexistence conditions, numerical solution strategy, parameter-estimation workflow, and software implementation. Section 3 demonstrates the framework through representative calculations, including idealized steady-state phase-diagram topologies, benchmark fitting to literature data, and reaction-dependent phase behavior in monodisperse and polydisperse systems. Additional derivations and implementation details are provided in the SI, including the chemical potentials, cloud-point equations and the Flory–Stockmayer distribution. The numerical results and example files used in the article are made available through the accompanying GitHub repository.⁴⁷ Section 4 summarizes the main conclusions, discusses the scope and limitations of the present implementation, and outlines directions for future work.

2. Methods

2.1 Overview of the PhaseTime framework

PhaseTime is a computational framework for predicting phase behavior and phase diagrams of quasibinary polydisperse polymer blends composed of two species (1 and 2) under protocol-dependent thermodynamic conditions. By combining Flory–Huggins thermodynamics^26,27 with modular representations of molecular weight distributions, interaction parameters, and experimental pathways, the framework enables the systematic computation of phase boundaries under both nonreactive and reactive conditions.

Fig. 1 illustrates the overall workflow of the PhaseTime framework. A central feature of the formulation is the representation of the thermodynamic state through a nondimensional scalar variable, denoted here as time t, which parametrizes the progression of the system along a prescribed experimental or computational protocol. If the protocol is explicitly time-dependent, t may represent physical time scaled by a reference value, t = t_phys/t_ref, or by a characteristic kinetic rate, for example t = kt_phys for a first-order process. More generally, t is used as a path variable that identifies the instantaneous state of the blend. This formulation is particularly useful for complex processing conditions, where the phase behavior is not controlled by temperature alone, but also by the evolving molecular constitution of the constituent species.

Fig. 1 PhaseTime framework: time-resolved evaluation of thermodynamic states and phase diagrams.

Within this framework, each state t is associated with a set of state functions. In the present implementation, these include the temperature,

T = T(t)

and extents of reaction of species 1 and 2,

a₁ = a₁(t), a₂ = a₂(t).

Isobaric conditions are assumed. It is further assumed that mass transport is fast relative to the chemical reactions, so that the system composition remains homogeneous, at least prior to phase separation. Accordingly, the thermodynamic state is taken to adjust instantaneously to the composition changes induced by the reaction process. Each one of the species 1 and 2 is considered as being derived from type-1 and type-2 monomers, respectively. Each monomer type may participate in chemical reactions with itself or with other constituents, forming new molecular constituents. However, no molecular constituents are derived jointly from both monomers 1 and 2. Typically, species 1 and 2 correspond to homopolymers of arbitrary architecture, formed by either step-reaction or chain polymerization (see sections 2.3 and 2.4). For given pressure and initial mixture composition, the quantities T(t), a₁(t), a₂(t) define the instantaneous thermodynamic and molecular state of the system; more generally, the state vector may be written as

s(t) = {T(t), a₁(t), a₂(t)}

from which all remaining model inputs are evaluated.

This protocol-based representation decouples the description of the experimental pathway from the underlying thermodynamic formalism. As a result, the same workflow can be applied to a wide range of scenarios, including isothermal reactive systems, nonreactive blends under temperature ramps, and coupled thermo-chemical protocols in which temperature and conversion evolve simultaneously. The state functions can be prescribed through arbitrary analytical expressions, thereby enabling flexible emulation of experimentally relevant schedules.

For each sampled state t, the framework evaluates the molecular distributions of products derived from the two species, the corresponding interaction parameter χ, and the thermodynamic conditions required for stability and phase equilibrium. In general, the interaction parameter may depend on temperature and/or composition through the protocol-dependent variables. Thus, the generalized parametrization by t provides a convenient means of tracing phase boundaries along arbitrary state trajectories without reformulating the thermodynamic problem for each specific protocol.

An important consequence of this formulation is that the resulting phase diagrams need not be restricted to the conventional temperature-composition representation. Since each point on a phase boundary is indexed by t, the vertical coordinate may be expressed in terms of any state-dependent quantity derived from the protocol, such as T, a₁, a₂, or χ. Likewise, composition may be represented using either φ₁ or φ₂. This flexibility allows the framework to generate phase diagrams in forms most appropriate for the system under study and the available experimental observables.

After prescribing a window of the generalized state variable t, at each point of the corresponding grid, the framework constructs the thermodynamic state of the system and evaluates the governing relations for phase stability and phase equilibrium. This enables the computation of spinodal curves, critical points, and coexistence boundaries, including the binodal in monodisperse blends and the cloud-point and shadow-point curves in polydisperse systems. The resulting sequence of states provides a unified description of phase behavior along the prescribed pathway, thereby enabling both prediction and direct comparison with experimental phase diagrams obtained under nontrivial processing conditions.

The framework supports fitting of model parameters to experimental phase-diagram data, including spinodal and coexistence curves, enabling both predictive modeling and quantitative interpretation of experimental phase behavior.

The framework is available both as a command-line workflow (CLI) and as a Python application programming interface (API). The command-line interface supports evaluation, fitting, residual diagnostics, plotting, and curve segmentation, while the Python API provides direct access to the underlying models and workflows through library imports. This dual design allows the framework to function both as an end-user tool and as a flexible library for custom scripting and method development.

2.2 Thermodynamic framework for quasibinary polydisperse blends and notation

Consider a quasibinary polymer blend comprising two polydisperse species, c ∈ {1, 2}. Each species c comprises constituents with n_ck molecules of degree of polymerization N_ck (counted in reference Flory–Huggins segments, each of volume V_ref), where k indexes the chain-size distribution for each species c; k ∈ {1, 2,…, k^(c)_max}. For notational convenience, the generic constituent index k is hereafter set to i for species 1 and j for species 2, following the notation adopted in ref. 6, 7 and 38. The total number of segments is:Assuming identical segment sizes, the volume fraction of the constituents c with chain size N_ck isThe volume fraction of each species is:Under the assumption that the system is incompressible,
Hereafter, the normalized chain-size (weight) fractions of species c will be denoted as w_ck and satisfy:The number-, weight-, and z-averages of species c, are defined with respect to the corresponding number of constituents or the corresponding volume fractions, as follows:Note that in the monodisperse limit, N_cn = N_cw = N_cz.

2.3 Flory–Huggins free energy and interaction models

2.3.1 Free energy and chemical potentials. In the context of the Flory–Huggins framework with temperature- and composition-dependent interactions,⁶ the normalized free energy of mixing per lattice site (also referred to as “free energy density” in the following) is:with Δ_mixF being the total free energy of mixing and g(T,φ₂) an interaction parameter dependent on composition and temperature, assuming a negligible dependence on chain size.³³ Given that the model is incompressible, with zero volume of mixing, Δ_mixF = Δ_mixA = Δ_mixG is both a Helmholtz energy and a Gibbs energy of mixing. The parameter g(T,φ₂) is commonly expressed in terms of an interaction parameter χ(T,φ₂) defined via chemical potentials;³⁸ the two are related as:^33,38

χ = g − g_φφ₁ (9)

where we have invoked the short-hand notation:

Upon integrating eqn (9),

The chemical potentials of the constituents 1i and 2j, relative to their pure state at the prevailing temperature and pressure, are defined as the partial derivatives of the free energy of mixing with respect to the number of molecules of these constituents, at constant temperature, pressure, and amounts of all other constituents:By substituting eqn (8) into eqn (11) and (12) we get:^6,7whereis the total dimensionless molecular-number density. Eqn (13) and (14) are identical to the expressions reported by Mumby et al.^6,7

We can reexpress eqn (13) and (14) in terms of χ(T,φ₂) as follows:

Analytical derivations of the chemical potentials can be found in section S1 (SI).

2.3.2 Models for the interaction parameter. We have implemented the functional forms based on ref. 38, where χ is a product of a temperature- and a composition-dependent term:

χ(T,φ₂) = D(T)B(φ₂) (18)

2.3.2/a Composition-dependent terms. According to ref. 38, the composition-dependent term may be represented adequately by a polynomial of the form

B(φ₂) = b₀ + b₁φ₂ + b₂φ₂² (19)

with b₀ = 1. The latter condition fixes the reference scale of the composition-dependent contribution and avoids introducing an unnecessary multiplicative degree of freedom. At the same time, the polynomial representation provides a convenient and flexible form for describing smooth variations of χ with composition. For this model, g(T,φ₂) in eqn (8) becomes:

which satisfies the relation in eqn (9).

2.3.2/b Temperature-dependent terms. Following ref. 38, the temperature-dependent factor D(T) is derived by separating it into enthalpic and entropic contributions.⁴⁸ Using the associated thermodynamic relations, and assuming a constant reduced heat capacity, integration yieldswhere T_ref = 1 K is introduced to render the logarithmic argument dimensionless. This expression provides a compact thermodynamically motivated representation of the temperature dependence of the interaction parameter. As shown in ref. 38 and related studies,^6,7 eqn (21) is sufficiently versatile to capture a broad range of phase-diagram topologies, including LCST, UCST, combined LCST-UCST behavior, as well as closed-loop and hourglass-type phase diagrams.

In addition, we consider the following polynomial representation:

for additional flexibility and testing purposes; e.g., setting D(T) ≡ T and T(t) = t, makes it possible to draw the phase diagrams in terms of a generalized interaction-strength variable D(T) as demonstrated in section 3.1.1.

2.4 Molecular distribution functions

The framework supports both monodisperse and polydisperse representations of the molecular constitution of the constituent species. These distribution models provide the set of chain sizes entering the thermodynamic formulation at each state t, and therefore establish the link between reaction progress, molecular structure, and phase behavior. In reactive systems, the distributions are evaluated from the prescribed extents of reaction, allowing the molecular populations to evolve along the prescribed protocol.

2.4.1 Monodisperse approximations. For applications in which a single effective chain size is sufficient, the framework supports monodisperse approximations where each species is represented by one characteristic degree of polymerization. In the present implementation, two such descriptions are implemented.

First, the framework implements a generalized Carothers-type expression in which the effective degree of polymerization is expressed as a function of the extent of reaction a_c of species c, according to:

This relation has been employed to represent the evolution of the number-average molar mass in the polycondensation of a stoichiometric diepoxy-diamine mixture upon setting γ_c = 4/3.^3,4 As such, it offers a simple analytical description of molecular growth and is especially useful in reduced reactive models, wherein the evolving molecular constitution is approximated by a single effective chain.

Second, the framework includes a monodisperse approximation of the Flory–Stockmayer description, where the full distribution is replaced by a single representative chain size identified with its analytically derived weight-average degree of polymerization, N_w. This approximation retains part of the conversion dependence implied by the underlying polymerization model, while avoiding the explicit treatment of the full chain-size distribution.

These monodisperse representations offer a computationally efficient route for coupling reaction progress to phase behavior, and are especially useful when the main objective is to capture the dominant shift of the phase boundaries with conversion rather than the full effect of polydispersity.

2.4.2 Polydisperse distributions. To account explicitly for conversion-dependent polydispersity, PhaseTime includes a generalized Flory–Stockmayer distribution model. In this formulation, a reactive species is represented as an initial mixture of multifunctional constituents carrying complementary A- and B-type reactive groups, which react exclusively with one another. The labels A and B therefore denote reactive functionalities within a given reactive species and should not be identified with species 1 and 2 of the quasibinary blend. The initial mixture may contain several classes of A-functional and B-functional constituents with different functionalities and reduced molecular volumes. As the reaction proceeds, the initial composition and the extent of reaction determine the evolving molecular population, yielding a conversion-dependent chain-size distribution.

For a given protocol state, the Flory–Stockmayer model^40,41 generates a discrete set of molecular constituents with reduced chain sizes N_k and normalized volume-fraction weights w_k. These quantities are then passed directly to the quasibinary Flory–Huggins formulation. Consequently, the evolving molecular distribution affects the entropy of mixing, the phase-stability conditions, and the coexistence boundaries. This representation therefore goes beyond an effective-chain approximation by retaining the explicit contribution of the molecular-size distribution at each state of the reaction protocol.

The moment evaluation follows the notation by Macosko & Miller⁴⁹ and Odle et al.,⁵⁰ in conjunction with the extensions proposed by Bachmann and Bendler,⁵¹ originally formulated for molecular-weight averages of cross-linked copolymers and branched polycondensates.^40,41 In PhaseTime, the same algebra is used after replacing molecular weights by reduced molecular volumes in Flory–Huggins lattice units. Thus, the resulting averages are interpreted as number-, weight-, and z-average reduced chain sizes, rather than true molar-mass averages.

The moments of the volume-weighted distribution are defined as

where k ∈ (1i, 2j), and the corresponding average reduced chain sizes are

The full Flory–Stockmayer distribution is infinite and must be represented numerically by a finite molecular population. The framework therefore constructs the distribution by enumerating molecular species up to prescribed precursor multiplicities and by discarding species whose relative volume fraction falls below a defined cutoff.

This finite representation is meaningful on the pre-gel side of the Flory–Stockmayer gel point. As the gel point is approached, the molecular-size distribution develops an increasingly extended long-chain tail and the analytical moments become singular; see sections S3.4 and S3.8 (SI). Consequently, increasingly large precursor multiplicities are required to represent the distribution accurately.

Since this truncation may remove a non-negligible contribution from the long-chain tail, an optional moment-correction procedure is provided. In this correction, effective tail species are introduced so that the analytical values of M₋₁, M₀, M₁, and M₂ are restored. This preserves N_cn, N_cw, and N_cz of the full distribution while retaining a finite and computationally tractable representation of the molecular population.

2.5 Phase stability and equilibrium boundaries

To facilitate the discussion and physical interpretation of phase boundaries under protocol-dependent conditions, we consider a representative reactive blend in which species 2 is monodisperse with degree of polymerization N₂ = 20, whereas species 1 evolves according to a Flory–Stockmayer chain-size distribution, initialized from A- and B-type precursors that react selectively with one another. The distribution was generated using A_fi = 0.4, B_gj = 0.6, N^A_i = 0.5, N^B_j = 0.8, f_i = 3, g_j = 2, m_max = 10, n_max = 10, with moment correction enabled; for a detailed description, please refer to section S3.3 (SI).

In the illustrative protocol considered here, temperature is prescribed as a function of nondimensionalized time as

and the extent of reaction of species 1 asThe first stage of the protocol corresponds to reaction under isothermal conditions, whereas the second combines a linear cooling ramp and further conversion with reduced rate. Finally, the temperature dependence of χ is described by eqn (21) with d₀ = 0.76, d₁ = −78 K and d₂ = 0 corresponding to LCST behavior; for simplicity, we treat χ as composition-independent; i.e., B(φ₂) = 1 in eqn (19).

Fig. 2c shows the prescribed evolution of T(t) and a₁(t), while Fig. 2d shows the corresponding evolution of the number-, weight-, and z-average degrees of polymerization, for species 1. The aforementioned quantities define the instantaneous state of the system at each value of t.

Fig. 2 Phase behavior under a prescribed protocol. Panels (a and b) show phase diagrams in the (φ₂, t) plane for the polydisperse system and its monodisperse counterpart, respectively. In both cases, the spinodal (black line) is identical, whereas the coexistence behavior differs: in the polydisperse system, distinct cloud and shadow curves are obtained (dashed and dot-dashed lines, respectively), which intersect the spinodal at the critical point (star); in the monodisperse limit, these curves collapse into a single binodal with α and β denoting its low-φ₂ and high-φ₂ branches, respectively. The cloud branches in panel (a) are denoted by α′ and β′, and the corresponding shadow branches by α″ and β″, respectively. Panel (c) shows the evolution of T(t) and a₁(t) and panel (d) the resulting number-, weight-, and z-average degrees of polymerization of species 1. The phase diagrams were computed and visualized using PhaseTime (based on Matplotlib). Additional curves were plotted using Veusz.⁵² Final figures were post-processed with Photopea.⁵³

To isolate the effect of polydispersity on phase behavior, we compare two systems. In the first, species 1 is described by the full Flory–Stockmayer distribution. In the second, species 1 is represented by a monodisperse approximation with degree of polymerization equal to the instantaneous weight-average value, N₁ = N_1w, dictated by the corresponding Flory–Stockmayer distribution. This comparison is particularly useful, because the spinodal depends only on the weight-average degrees of polymerization and is therefore identical in the two cases, whereas the coexistence curves differ due to the presence/absence of polydispersity.

Fig. 2a and b present the resulting phase diagrams in the (φ₂, t) plane. Although the vertical axis is expressed in terms of t, each point on the diagram corresponds to a specific thermodynamic state defined by the imposed temperature and reaction protocol. Both diagrams exhibit stable, metastable, and unstable regions, indicated by the different shades. The boundary between the unstable and metastable regions is defined by the spinodal (thick black line), while the boundary between the metastable and stable regions is determined by the coexistence conditions.

In the polydisperse system, the onset of phase separation is described by the cloud-point curve (dashed line), while the corresponding shadow-point curve (dot-dashed line) gives the composition of the infinitesimal coexisting phase that first appears at the cloud point. These curves meet the spinodal at the critical point (star). In contrast, for the monodisperse reference system, the cloud and shadow curves collapse into a single curve, the binodal. The side-by-side comparison therefore highlights directly the effect of polydispersity on phase equilibrium, while preserving the same underlying spinodal boundary.

Hereafter, coexistence quantities are grouped into an α-family and a β-family, corresponding to the low-φ₂ and high-φ₂ sides of the diagram, respectively. For monodisperse systems, the binodal branches are denoted simply by α and β, accordingly; quantities associated with these branches are identified by the superscripts (α) and (β), respectively. For polydisperse systems, quantities associated with the cloud phase are denoted by a prime superscript (x′), whereas those associated with the shadow phase are denoted by a double-prime superscript (x″), following the notation in ref. 7. Thus, α′ and β′ in Fig. 2a denote the cloud branches, while α″ and β″ denote the corresponding shadow branches.

2.5.1 Spinodal curve. The spinodal is obtained from the condition that the Hessian of the free energy of mixing becomes singular; i.e., the determinant of the Hessian (H_mix) with respect to the composition variables vanishes,
where ij and kl run over all constituents of the two species; conveniently, this condition can be expressed explicitly as

In the present comparison, the spinodal is identical for the polydisperse system and for its monodisperse N₁ = N_1w counterpart, since the spinodal depends only on the weight-average degrees of polymerization of the two species (cf. Fig. 2a and b).

2.5.2 Critical points. The critical state is determined by simultaneously satisfying the spinodal condition eqn (28) and the corresponding higher-order criticality criterion:

Y′ = 0 (29)

with Y′ denoting the determinant obtained from eqn (28) by replacing the elements of any one row with the corresponding third derivatives.^6,33 Eqn (29) can be expressed explicitly as follows:When these conditions are solved for each prescribed state t, they define a critical line in the (φ₂, t) plane, i.e., the locus of critical compositions associated with the instantaneous thermodynamic state and molecular distribution imposed by the protocol.

The critical point observed in the phase diagram is therefore a particular point on this critical line, corresponding to the state at which the coexistence boundary intersects the spinodal. The full critical line, however, provides a more complete description, since it reveals how the critical composition evolves as the molecular constitution of species 1 changes during the protocol. In the present case, this evolution reflects the coupled effects of the imposed thermal pathway and the conversion-dependent changes in the Flory–Stockmayer distribution.

Comparison of the polydisperse system with its monodisperse reference counterpart further clarifies the role of polydispersity. Although the two systems share the same spinodal, because this depends on the weight-average degrees of polymerization, their critical lines do not in general coincide; cf. Fig. 2a and b. The difference arises from the coexistence conditions, which are altered by the full molecular distribution in the polydisperse case. The critical lines shown in the phase diagrams therefore provide a direct visualization of how polydispersity shifts the location of criticality, beyond its effect on the coexistence branches themselves. In the monodisperse limit, however, the critical point coincides with the extremum of the spinodal curve, reflecting the fact that the coexistence and stability boundaries become tangential at that point.

Note that, for the special case in which the interaction parameter is composition-independent, the extremum and the critical composition of a quasibinary polydisperse blend admit the following analytical expressions,

whereas the corresponding χ_extr and χ_crit are obtained by substituting φ_2,extr and φ_2,crit, respectively, into the spinodal condition, eqn (28). In the monodisperse limit where N_cw = N_cz = N_c, this reduces to the classical Flory–Huggins result, φ_2,crit = φ_2,extr. Although the unique spinodal extremum and the critical point coincide in the classical monodisperse case, they are treated as distinct objects in the general quasibinary polydisperse framework.

2.5.3 Binodal curves. The coexistence boundary marks the onset of phase separation, separating the single-phase region from the two-phase region. In monodisperse blends, this boundary is defined by the binodal curve, which corresponds to the locus of points in the (φ₂, t) plane where the chemical potentials of the two species are equal, i.e., Δμ^(α)₁ = Δμ^(β)₁ and Δμ^(α)₂ = Δμ^(β)₂, where eqn (16) and (17) give:for the present model.³⁸

2.5.4 Cloud- and shadow-point curves. In polydisperse systems, the coexistence boundary is described by cloud- and shadow-point curves. At the cloud point, a prescribed homogeneous parent phase first coexists with an infinitesimal conjugate phase. The parent phase defines the cloud point, whereas the infinitesimal conjugate phase defines the corresponding shadow point. In PhaseTime, these curves are evaluated using the approach of Mumby et al.,⁶ which extends an earlier implementation⁵ based on the treatments of Koningsveld and co-workers³³ and Šolc³⁴ from quasi-binary systems with one polydisperse component to systems in which both components are polydisperse.⁶ Their equilibrium is determined by equality of the chemical potentials of all constituents in the two phases:In the general case, the coexistence problem requires solving a coupled nonlinear system whose dimensionality grows with the number of molecular constituents used to represent the two polydisperse species. However, when the relative molar-volume distributions of the two polymers are known prior to phase separation, the problem can be reduced substantially by introducing two separation factors, σ₁ and σ₂. Within the Flory–Huggins framework, these factors describe the partitioning of each constituent between the cloud and shadow phases and provide a compact representation of fractionation during phase separation:

At the cloud point, one of the coexisting phases coincides with the pre-existing homogeneous parent phase, here termed the principal phase, whose molecular distribution is known. The composition of the second, infinitesimal shadow phase can then be obtained from the partitioning relations. As a result, the full many-constituent coexistence problem collapses to a reduced system involving the state variable, the total composition of the principal phase, and the separation factors. In particular, substituting eqn (37) and (38) into eqn (35) and (36) yields the following system of equations:

We note that the sign in eqn (39) is opposite to that reported in ref. 6 and 7; this discrepancy may arise from a difference in notation or a typographical error in the original expressions. An analytical derivation of eqn (39) and (40) is provided in section S2 (SI).

In the present protocol-dependent formulation, these coexistence conditions are evaluated at each prescribed state t, with the molecular distribution of species 1 determined by the corresponding reaction extent. The cloud and shadow branches shown in Fig. 2a therefore represent the onset of demixing and the composition of the associated incipient phase, respectively, along the imposed thermo-reactive pathway. Their intersection with the spinodal identifies the critical point, while comparison with the monodisperse reference system in Fig. 2b illustrates that, in the absence of polydispersity, the cloud and shadow curves collapse into a single binodal curve. This comparison highlights directly the role of polydispersity in broadening the coexistence description from a single binodal into distinct cloud and shadow branches.

2.6 Numerical workflow and solution strategy

2.6.1 Protocol sampling and state construction. For a prescribed protocol, the framework treats the generalized state variable t as the independent sampling coordinate and constructs a discrete grid {t_k} over a prescribed interval. At each grid point, the framework evaluates the corresponding thermodynamic state, including the instantaneous temperature, protocol variables, molecular descriptors of the two quasibinary species, and the interaction parameter χ. This design decouples the numerical solver from the specific form of the experimental or reactive protocol, allowing arbitrary user-defined state trajectories to be treated within the same computational workflow.

All phase-boundary calculations are performed from these protocol-defined states. Depending on the requested outputs, PhaseTime evaluates the spinodal, critical features, binodal (monodisperse approximation), and cloud/shadow branches in a modular manner, and subsequently assembles the resulting data into phase-boundary datasets.

2.6.2 Spinodal evaluation. For each sampled value of the generalized state variable t, PhaseTime constructs the instantaneous thermodynamic state and evaluates the corresponding spinodal function χ_sp(T,φ₂), eqn (28).

Thereafter follows the determination of the stationary points of this function, i.e. the compositions at which

These stationary points define the extrema of the spinodal function for the instantaneous state, whereas the lowest extremum defines the minimum value of χ_sp required for the existence of real solutions of eqn (28). The extrema are determined by evaluating eqn (41) using a model-specific interaction routine that accounts for the composition dependence of the interaction parameter; i.e. B(φ₂) in eqn (18).

In case the protocol-defined interaction parameter satisfies

χ(t) ≥ χ_sp,min(t),

the spinodal branches are obtained by solving eqn (28) over 0 < φ₂ < 1, by incorporating Brent's method⁵⁴ as implemented in SciPy.⁵⁵ This formulation is general and accommodates models for which the spinodal function may possess more than one stationary point.

2.6.3 Detection of Extrema and critical points. A refined determination of the spinodal extrema and critical points is performed in a separate stage after the state-wise spinodal sweep. The extrema are evaluated for the instantaneous state from the conditions of eqn (41). The lowest extremum, which defines the onset of the two-root spinodal regime along the protocol, is then identified by solving the equation,

χ_sp,min(t) − χ(t) = 0 (42)

Critical points are evaluated independently as follows:

χ_crit(t) − χ(t) = 0 (43)

where χ_crit is obtained by solving eqn (28) and (30). The solution strategy for eqn (42) and (43) depends on whether χ is composition-dependent; if χ is composition-independent, analytical solutions may be obtained [i.e., see eqn (31) and (32)], otherwise the equations are solved numerically.⁵⁵

2.6.4 Binodal continuation. The binodal is computed by solving the coexistence conditions for a pair of coexisting phase compositions and then continuing the solution along the generalized state variable t. At each protocol state, the unknown phase compositions are obtained by enforcing equality of the two quasibinary chemical potentials between the coexisting phases in eqn (33) and (34), or equivalently by finding pairs of (φ^(α)₂, φ^(β)₂) for which the following residuals in eqn (44) and (45) vanish.

r₁(φ^(α)₂,φ^(β)₂) = Δμ^(α)₁ − Δμ^(β)₁ (44)

r₂(φ^(α)₂,φ^(β)₂) = Δμ^(α)₂ − Δμ^(β)₂ (45)

The nonlinear system is solved by bounded least-squares optimization using the Trust Region Reflective (TRF) algorithm in SciPy,⁵⁵ following the bound-constrained trust-region formulation of Branch, Coleman, and Li.⁵⁶

Initial coexistence pairs may be generated in three ways: i) by scanning a two-dimensional composition grid, ii) by applying prescribed offsets to the local spinodal roots, or iii) by direct specification. In the scan-based approach, the search domain can be substantially reduced without loss of distinct coexistence solutions by excluding the trivial diagonal and the spinodally unstable region (see Fig. 3).

Fig. 3 Residual landscape of the binodal eqn (44) and (45), F(φ^(α)₂,φ^(β)₂) = 0.5(r₁² + r₂²), evaluated across the (φ^(α)₂,φ^(β)₂) grid at a representative protocol state (t = 60, Fig. 2b). Thick solid lines indicate the spinodal compositions. The symmetry about the diagonal φ^(α)₂ = φ^(β)₂ reflects the trivial solutions, while the dashed green rectangle marks the reduced domain used for initialization of the binodal search, excluding both symmetric duplicates and the spinodally unstable region. The red circles mark the non-trivial zeros for (φ^(α)₂,φ^(β)₂) = {(0.13,0.32), (0.32,0.13)}.

Each accepted seed is then continued forward and/or backward in t by using the previously converged coexistence pair as the initial guess at the next protocol state. If convergence fails, the continuation step is reduced adaptively until the solution is recovered or a minimum step size is reached. Starting states may also be introduced automatically from intermediate locations between special points (extrema and critical spinodal points and protocol boundaries), thereby improving branch detection over the full protocol range. The final binodal dataset is assembled from all successful continuation traces.

2.6.5 Cloud- and shadow-point continuation. The cloud and shadow curves are computed by a continuation procedure analogous to that used for the binodal, but specialized for incipient phase separation in the polydisperse system. At a given protocol state, the unknowns are the cloud-point composition of the parent phase, , and an auxiliary reweighting parameter, σ₂, from which the infinitesimal shadow phase is reconstructed. Here, and denote the normalized chain-size fractions of species 1 and 2 in the parent phase, whereas and denote the corresponding fractions in the shadow phase; see eqn (4).

For species 2, the shadow molecular-weight distribution is obtained by exponential reweighting of the parent distribution,

so that the corresponding shadow composition is
An analogous relation is used for species 1 through a second reweighting parameter, σ₁.

For a trial value of , the shadow phase is reconstructed self-consistently by first evaluating , then determining σ₁ from the first cloud eqn (39), and finally obtaining the shadow distribution and composition of species 1. In detail, σ₁ is eliminated from eqn (39) as follows:

after which
thus,
The cloud solution is then obtained by enforcing the remaining two conditions: the incompressibility constraint,and satisfying the second cloud eqn (40),therefore, the cloud-point problem reduces to the determination of for which both residuals vanish. Numerically, the nonlinear system is solved by bounded least-squares optimization.^55,56 Accepted solutions are required to satisfy successful evaluation of the cloud equations, numerical convergence of the optimizer, finite separation between cloud and shadow compositions in order to exclude the trivial solution, and equality of the constituent chemical potentials between parent and shadow phases within a prescribed tolerance.

As in the binodal calculation, the procedure is initialized from one or more starting points at selected protocol states. These starting points may be obtained either from a scan over or from direct user specification. In the scan-based initialization, candidate solutions are sought separately on the low- and high-composition sides of the critical composition, and only admissible converged solutions are retained; e.g., see Fig. 4. Each accepted seed is then continued forward and/or backward along the generalized state variable t by reusing the previously converged pair as the initial guess at the next state. If convergence fails, the continuation step is reduced adaptively until the solution is recovered or a minimum step size is reached. At each successful step, the cloud-point composition of the parent phase and the corresponding shadow-point composition of the infinitesimal phase are recorded separately, yielding paired cloud and shadow branches along the protocol.

Fig. 4 Residual landscape, from eqn (47) and (48) on grid at a representative protocol state (t = 60, Fig. 2a). Thick solid lines indicate the spinodal compositions. The line σ₂ = 0 corresponds to the trivial solution where , while the dashed green rectangles mark the optimal reduced domains for scan-based initializations after exclusion of the spinodally unstable region. Red circles denote the non-trivial zeros of the cloud residuals for this protocol state, .

2.7 Scope and limitations of the thermodynamic framework

Before discussing parameter estimation and practical use of the software, we summarize the main assumptions and limitations that define the scope of the present thermodynamic model.

2.7.1 Assumptions of mean-field lattice thermodynamics. As a Flory–Huggins-based framework, PhaseTime inherits the assumptions of mean-field lattice thermodynamics. Its main limitations concern: (i) the form of the free energy, which does not explicitly resolve conformational or packing effects; (ii) the molecular populations that can be represented; and (iii) the use of a lattice-based volume mapping and incompressibility assumption.

The present model does not explicitly describe polymer stiffness, chain flexibility, persistence length, orientational ordering, or conformation-dependent entropy. Stiffness-induced effects may be incorporated only phenomenologically, for example through an effective Flory–Huggins interaction parameter as proposed for stiffness-mismatched blends by Kozuch et al.⁵⁷ However, the present formulation cannot predict phase separation driven directly by conformational or orientational degrees of freedom, such as the nematic unmixing of semiflexible polymers reported by Milchev et al.⁵⁸ Such effects would require an extended free-energy model with explicit conformational, orientational, or architecture-dependent contributions.

The current implementation is restricted to quasibinary blends in which the two polymeric species remain distinct. Molecular growth may occur within each species, but the framework does not generate copolymer molecules containing units from both species 1 and 2. A fixed-composition copolymer could be treated phenomenologically as an effective component with fitted parameters, but this would not generally be transferable to other compositions or architectures. Explicit random, block, graft, or reactive copolymer formation would require both an extended distribution model and a corresponding free-energy description that tracks copolymer composition and architecture.⁵⁹ Furthermore, the implemented distribution models do not explicitly account for intramolecular cyclization or cyclic polymer species. In particular, the Flory–Stockmayer distribution assumes tree-like branching^40,41 and is appropriate only when ring formation is negligible or can be absorbed into effective distribution parameters. More refined distribution models could in principle be incorporated into the framework,^60,61 but cyclization-corrected distributions are not included in the present implementation.

The use of a single reference volume is part of the coarse-grained lattice approximation. The reference volume V_ref defines reduced chain sizes, N_i = V_i/V_ref. Its numerical choice is not unique; changes in V_ref can be absorbed by the corresponding rescaling or refitting of the interaction parameter. Thus, blends with different molecular or segment sizes can be treated phenomenologically. However, in doing so, size and packing effects are represented through reduced chain sizes and effective interaction parameters, rather than through an explicit microscopic packing description.

Finally, the present formulation is built upon the incompressibility constraint (φ₁ + φ₂ = 1) which is appropriate for dense polymer blends and melts where volume changes upon mixing and pressure-volume effects are small. Deviations from incompressibility, such as density changes, free-volume effects, or pressure-dependent miscibility, may be treated phenomenologically to some extent through fitted effective interaction parameters. Nevertheless, the model does not include an explicit density or free-volume degree of freedom; systems in which compressibility or volume changes are central would therefore require a compressible extension, for example through an equation-of-state description.

2.7.2 Quasi-equilibrium interpretation of phase diagrams. The quasi-equilibrium interpretation of the protocol-dependent phase diagrams requires composition relaxation to be fast compared with the reaction or imposed protocol. A useful order-of-magnitude estimate is obtained by comparing the demixing time, τ_demix ∼ l²/D_rel with a characteristic protocol timescale, τ_protocol. Here, l is a representative composition-relaxation length scale, such as the dominant spinodal wavelength,⁶² droplet size, or domain size. D_rel is an effective mutual diffusion coefficient; depending on the kinetic description, it may be estimated from a mutual/interdiffusion coefficient⁶³ or from the mobility function employed in phase-field models.⁶⁴ The characteristic protocol time τ_protocol may correspond to a reaction time scale, for example τ_rxn ∼ Δa/|da/dt|, or to a temperature-ramp time scale, τ_T ∼ ΔT/|dT/dt|.

For representative sub-micron/micron morphologies, l = 0.1–1 μm,⁶⁵ and polymeric diffusion coefficients D_rel = 10⁻¹⁴–10⁻¹² m² s⁻¹, τ_demix is approximately 10⁻²–10² s. For slower diffusion, D_rel = 10⁻¹⁶–10⁻¹⁵ m² s⁻¹, relaxation over the same length scales can range from seconds to hours. Thus, the quasi-equilibrium approximation is system- and protocol-dependent and may fail near vitrification, gelation, or late-stage curing. In such cases, the present calculations should be used as a thermodynamic reference for kinetic reaction–diffusion treatments and for phase-field models describing steady-state and reaction-induced phase separation (RIPS).

2.8 Parameter estimation from experimental phase-diagram data

PhaseTime provides a fitting workflow for estimating interaction-model parameters from experimental phase-diagram data. Experimental input is provided as discrete (φ₂, t) points sampled along a phase boundary, with the generalized state variable t defined, as in the general evaluation workflow, to represent the evolution of temperature (T), extents of reaction (a₁, a₂), or mixed protocol conditions (see section 2.1). At present, the fitting module supports spinodal datasets and binodal-type coexistence datasets represented through the formulation in section 2.5. For each trial parameter vector, the framework reconstructs the model and reevaluates the corresponding phase boundary using the numerical machinery described in section 2.6.

2.8.1 Objective function. The experimental curve is preprocessed once prior to optimization and, where required, decomposed into N^exp_b distinct single-valued t(φ₂) segments, to enable comparison of multibranch phase boundaries. Details of the segmentation algorithm are given in section S4 (SI).

For each trial parameter vector θ, the model is reconstructed and the target phase boundary is reevaluated over the prescribed grid of the generalized state variable t, augmented with the experimental t-values. The computed phase boundary is then decomposed into N^mod_b single-valued branches and compared with the experimental dataset. For each matched branch pair, the comparison is restricted to their common composition interval, and both branches, represented as t(φ₂), are linearly interpolated onto a shared set of φ₂-evaluation points. The branchwise mismatch is then defined as the mean squared normalized residual in t, while an additional penalty is introduced when the computed branch does not cover the full experimental composition span. Solutions that yield a different number of branches than the experimental dataset, or otherwise fail during evaluation, are assigned a large penalty value Λ.

In the present implementation, we introduce a normalization factor s_t based on a global characteristic scale of the experimental dataset, taken as the overall span of the generalized state variable t across the full dataset:

s_t = max t^exp − min t^exp

so that the residuals are scaled by the overall extent of the experimental curve in t.

The corresponding objective function may be written as

The branchwise shape mismatch is defined as

where N_k is the number of φ₂-evaluation points for branch k, whereas t̃^exp_k,i and t̃^mod_k,i are the experimental and computed values obtained by linear interpolation onto the common composition grid. The term
penalizes incomplete coverage of the experimental branch by the computed branch in composition space, increasing quadratically with the missing fraction of the experimental φ₂-span. Here, Δφ^exp_k is the composition span of the experimental branch, Δφ^overlap_k is the overlap span between experimental and computed branches, and w_φ is a fixed penalty coefficient. In this way, the objective function rewards agreement in branch shape and extent, while strongly penalizing topologically inconsistent or numerically invalid solutions.

2.8.2 Optimization workflow. The current implementation supports both local and derivative-free optimization. Local optimization is performed through SciPy minimization routines⁵⁵ with bound handling, whereas derivative-free search may be carried out using the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) algorithm.^66,67 In both cases, the workflow returns the optimal parameter vector, the final loss, the updated model configuration, and the optimization history. The present parameter-estimation module should be regarded as a preliminary constituent of PhaseTime, intended primarily for exploratory fitting studies. Nevertheless, the current implementation has already been tested on the fitting of several spinodal and binodal curves exhibiting LCST, UCST, LCST + UCST, hourglass, and loop topologies. Future development will focus on improving robustness and convergence reliability across a broader range of datasets, including cloud- and shadow-point curves.

2.9 Software implementation and user workflow

PhaseTime is a Python-based computational framework that supports both end-user operation through a command-line interface and direct programmatic use via a Python API. Both interfaces share a common computational backend and configuration logic, ensuring consistency and reproducibility across workflows. This dual design enables routine tasks—such as phase-diagram evaluation, fitting, plotting, and diagnostic analysis—while also supporting more customized computational workflows for method development. The following sections summarize the modular software architecture underlying the framework and the main user interfaces.

2.9.1 Software architecture and extensibility. PhaseTime is structured as a modular framework in which thermodynamic models, molecular-distribution models, protocol functions, and workflow orchestration are implemented as separate constituents. A central factory routine assembles these constituents from the specified configuration into a unified model object, allowing the same backend to support both the command-line interface and the Python API.

The composition- and temperature-dependent parts of the interaction parameter are implemented through separate abstract model hierarchies with registry-based construction. New model forms can therefore be introduced by defining and registering an additional subclass, without changing the higher-level evaluation workflow. Molecular distributions are handled analogously through dedicated distribution functions operating on a common polymer representation with cached molecular-weight averages.

This organization promotes extensibility, code reuse, and consistency across evaluation, fitting, and diagnostic tasks, and facilitates the incorporation of additional constitutive relations or numerical procedures as the framework evolves.

2.9.2 Command-line workflow. PhaseTime provides a configuration-driven command-line interface for model evaluation, parameter estimation, visualization, and numerical diagnostics. The main executable (phase-time) exposes dedicated subcommands for phase-diagram evaluation (eval), fitting to experimental data (fit), and plotting of previously computed results (plot). Auxiliary utilities are also provided for diagnostic analysis. The subcommand binodal-grid evaluates the residuals of the binodal eqn (44) and (45) over a prescribed grid, while the subcommand cloud-grid evaluates the residuals of the cloud-point eqn (47) and (48) in an analogous manner; these tools enable visualization of the residual landscapes, as illustrated in Fig. 3 and 4, respectively. Finally, the subcommand segment-curve preprocesses input phase-diagram data by decomposing multivalued curves into one or more single-valued segments, enabling consistent visualization.

In routine use, PhaseTime is configured through an external YAML input file that specifies the thermodynamic model, molecular distributions, protocol functions, state grid, and numerical settings. This configuration-driven structure separates model specification from execution and supports reproducible calculations.

2.9.3 Python API. PhaseTime provides a Python API that exposes the same underlying thermodynamic and numerical functionality for programmatic use. This interface enables direct access to phase-diagram evaluation, fitting, and diagnostic routines from Python scripts or interactive environments, thereby supporting custom workflows beyond the standard command-line use cases. The Python API is intended for applications such as automated parameter studies, integration into larger simulation or optimization pipelines, custom post-processing, and user-defined visualization. By returning the computed phase-boundary data directly to the Python environment, it facilitates flexible downstream analysis without requiring intermediate manual steps.

3. Results

This section demonstrates the capability of PhaseTime to reproduce and analyze phase behavior in quasibinary polymer blends under both quiescent and reaction-dependent conditions. The results are organized to first assess equilibrium phase diagrams in nonreactive systems, where the framework is tested against representative monodisperse and polydisperse cases, and then to examine protocol-dependent phase behavior during reaction. Emphasis is placed on the ability of the framework to capture diverse coexistence-topology classes, incorporate evolving molecular distributions, and resolve the effect of distribution truncation and reduced-order approximations on the predicted phase behavior.

3.1 Steady-state phase behavior

3.1.1 Idealized phase-diagram topologies in the monodisperse limit. We first consider monodisperse blends under steady-state conditions to illustrate the range of phase-diagram topologies accessible within the present framework. For validation, representative UCST, LCST, combined UCST-LCST, hourglass, and loop diagrams were evaluated using the models defined by eqn (19) and (21) in conjunction with the parameter sets listed in Table 1, from ref. 38.

Table 1 Parameters of the temperature- and composition-dependent contributions to the interaction parameter, χ, defined by eqn (21) and (19), respectively, for the monodisperse blends with N₁ = 100 and N₂ = 10 000. Units: d₁/K

Case	b1	b2	d0	d1	d2	Type
a	0	0	0.013	−2.01	0	LCST
b	0.01	0.001	0.013	−2.01	0	LCST
c	0.01	0.001	0	2.01	0	UCST
d	0.01	0.001	−0.03487	2.01	0.006	UCST + LCST
e	0.01	0.001	−0.03435	2.01	0.006	Hourglass
f	0.01	0.001	0.047	−2.01	−0.006	Closed-loop

---

The topology is governed mainly by the signs of d₁ and d₂.³⁸ For d₂ = 0, d₁ > 0 gives UCST behaviour, whereas d₁ < 0 gives LCST behaviour. When d₂ ≠ 0, the additional temperature dependence permits nonmonotonic coexistence: d₁ > 0, d₂ > 0 yields combined UCST-LCST or hourglass topologies, whereas d₁ < 0, d₂ < 0 produces a closed loop. Thus, d₁ sets the primary temperature trend of the interaction parameter, and d₂ controls its curvature. The calculated phase diagrams are shown in Fig. 5 and follow these expected trends.

Fig. 5 Calculated phase diagrams for monodisperse blends with N₁ = 100 and N₂ = 10 000. Panel labels correspond to the cases listed in Table 1.

It should be noted that the coefficients controlling the strength of the composition-dependent term, b₁ and b₂, are small in Table 1. This is illustrated by the case in Fig. 5a with no composition dependence, which is very similar to its composition-dependent counterpart in Fig. 5b. Therefore, in the following, we will showcase an example with interaction parameters which are strongly composition-dependent, leading to exotic behavior which can be found in quasibinary systems.

Depending on the functional form of the interaction parameter, the critical-point equations may admit multiple solutions. This possibility is illustrated in Fig. 6 for a symmetric monodisperse blend with N₁ = N₂ = 10 and

χ(t,φ₂) = tB(φ₂), B(φ₂) = 1 − 3φ₂ + 2.8φ₂² (50)

where t ≡ D(T) is used as a generalized interaction-strength variable. Solving the critical-point eqn (28) and (30) gives three critical compositions: two local minima of the spinodal, (t_c, φ_2,c) = (0.31, 0.82) and (0.87, 0.14), and one local maximum, (t_c, φ_2,c) = (2.17, 0.39).

Fig. 6 Phase behavior of a symmetric monodisperse blend with N₁ = N₂ = 10 and interaction parameter from eqn (50). The prefactor t ≡ D(T) is used as a generalized interaction-strength variable. The spinodal exhibits three critical points, corresponding to two local minima and one local maximum. Multiple common-tangent branches are obtained from the coexistence equations; the globally stable binodal is selected by the lower-convex-envelope criterion, while the additional branch corresponds to hidden/local coexistence.

The coexistence equations consequently admit multiple common-tangent branches. However, these branches are not all thermodynamically equivalent. The physically stable binodal is the branch whose tie lines form supporting lines of the free-energy density, or equivalently the branch selected by the lower convex envelope of f_m(φ₂; t). Additional branches may satisfy the local common-tangent equations, but are preempted by the globally stable two-phase state and should therefore be interpreted as hidden or local coexistence branches rather than equilibrium binodals.

This distinction is clarified in Fig. 7, which shows the corresponding free-energy density f_m(φ₂; t) at selected values of t, together with the relevant tangent constructions. At the lowest critical point, t = 0.31, the globally stable tie line has collapsed to a single tangent at φ_2,c = 0.82. For a nearby value inside the coexistence region, e.g. t = 0.50, a finite common tangent connects two distinct coexisting compositions and forms a supporting line of the free-energy curve. This is the usual equilibrium binodal construction.

Fig. 7 Free-energy density f_m(φ₂, t) from eqn (8) for the model shown in Fig. 6. Orange curves show f_m(φ₂, t), and blue lines show the corresponding tangent constructions. Vertical dashed/dotted lines mark the compositions of global/hidden critical or coexistence points. At t = 0.31, the stable tie line collapses to the lower critical point, while at t = 0.50 a finite supporting common tangent defines the stable binodal. At t = 0.87 and t = 2.17, the dashed blue lines show collapsed tangents associated with hidden critical points, whereas the solid blue lines indicate the globally stable supporting tangents. Only the latter form the lower convex envelope and correspond to the equilibrium binodal.

At t = 0.87, the local inner branch reaches a critical point at φ_2,c = 0.14, so its tie line collapses to a tangent at that composition. However, the same free-energy curve also admits a lower outer supporting tangent connecting compositions close to the composition boundaries. Therefore, the local critical point is hidden by the globally stable two-phase state. The same interpretation applies at t = 2.17, where the upper local critical point at φ_2,c = 0.39 corresponds to a collapsed tangent of the hidden branch, while the globally stable coexistence is still determined by the outer supporting line.

Thus, the presence of three critical points does not imply three independent globally stable binodals. Rather, the coexistence equations define several local common-tangent branches, and the physical binodal is obtained only after applying the convex-envelope criterion. This provides a direct numerical validation that the additional branches are genuine roots of the coexistence equations, but not all of them correspond to equilibrium phase boundaries.

3.1.2 Literature benchmarks for monodisperse phase-boundary fitting. The monodisperse framework was fitted to the experimental phase-boundary data listed in Table 2. The selected systems span several classes of phase behaviour and include both binodal and spinodal datasets. The phase data^68–72 were taken from the literature datasets compiled in ref. 38 with the original experimental sources cited in Table 2 for each case. These cases are used here as benchmark problems for assessing the numerical implementation and flexibility of PhaseTime, rather than as a reanalysis of the underlying experiments.

Table 2 Molecular weights, mass densities, chain molar volumes, and dataset types for the experimental systems considered in the monodisperse analysis. Subscripts 1 and 2 denote the first and second components in the case label, respectively. Trailing suffixes indicate molecular-weight labels or distinct system variants. The abbreviations bn and sp denote binodal and spinodal datasets, respectively. Units: M_c/g mol^–1, ρ_c/g cm^–3, V_c/cm³ mol^–1

Case	M1	ρ1	V1	M2	ρ2	V2	Description
PVME_PS62k	75 000	1.000	75 000	62 000	1.000	62 000	LCST (bn)^69
PVME_PS100k	75 000	1.000	75 000	100 000	1.000	100 000	LCST (bn)^69
PVME_PSD_L	389 000	0.983 (ref. 38 and 73)	395 727	230 000	1.092 (ref. 38 and 74)	210 623	LCST (sp)^68
PVME_PSD_H	593 000	0.983 (ref. 38 and 73)	603 255	1 100 000	1.092 (ref. 38 and 74)	1 007 326	LCST (sp)^68
PBD_DPBD	135 000	0.884 (ref. 70)	152 715	134 000	0.982 (ref. 70)	136 456	UCST (sp)^70
acetone_PS5k	58	0.7845 (ref. 75)	74	4800	1.060 (ref. 76)	4528	LCST + UCST (bn)^71
acetone_PS20k	58	0.7845 (ref. 75)	74	19 800	1.060 (ref. 76)	18 679	Hourglass (bn)^71
water_PVA	18	1.000	18	140 000	1.300 (ref. 77)	107 692	Loop (bn)^72

---

The digitized⁷⁸ phase-boundary data were refitted using an effective monodisperse description of each component. Molecular weights, mass densities, and related system properties were taken from the original articles and references therein, where available. However, the reference volume V_ref used to define the Flory–Huggins segment scale was not reported in ref. 38. Therefore, the chain sizes were recalculated using the reference volumes listed in Table 2, according to

with M_c and ρ_c being the effective molecular weight and mass density of species c, respectively. Since V_ref sets the scale of the segment-based chain sizes and interaction parameters, the model parameters were refitted within PhaseTime rather than copied directly from ref. 38. Consequently, some fitted coefficients differ slightly from those reported in ref. 38, while yielding the same phase-boundaries within the adopted parameterization.

For systems where specific density values or reference volumes could not be identified, representative density values were assigned from the literature; in the absence of more specific information, a density of 1 g cm⁻³ was used. The reference volume affects the numerical scale of the fitted interaction parameters but not the phase boundary itself, provided the fitted coefficients are rescaled consistently. The adopted reference volumes, resulting chain sizes, and fitted interaction-parameter coefficients are reported in Table 3.

Table 3 Flory–Huggins reference volumes, effective chain sizes, and fitted interaction-parameter coefficients used in the monodisperse calculations for the systems listed in Table 2. The chain sizes N₁ and N₂ are dimensionless and were computed according to eqn (51). Units: V_ref/cm³ mol^–1, d₁/K, everything else is dimensionless

System	Vref	N1	N2	b1	b2	d0	d1	d2
PVME_PS62k	58.08	1291	1067	−3.640 × 10^−1	0	1.709 × 10^−2	−6.000	0
PVME_PS100k	58.08	1291	1722	−3.640 × 10^−1	0	1.709 × 10^−2	−6.000	0
PVME_PSD_L	59.08	6698	3565	−1.605	8.444 × 10^−1	4.098 × 10^−2	−1.673 × 10	0
PVME_PSD_H	59.08	10 210	17 049	−1.706	9.851 × 10^−1	2.027 × 10^−2	−8.433	0
PBD_DPBD	61.10	2499	2233	−5.284 × 10^−1	2.114 × 10^−1	8.988 × 10^−5	2.480 × 10^−1	0
acetone_PS5k	74.03	1	61	2.031 × 10^−1	0	−6.799	3.524 × 10^2	1.083
acetone_PS20k	74.03	1	252	6.962 × 10^−1	0	−7.009	3.753 × 10^2	1.103
water_PVA	18.02	1	6000	6.542 × 10^−1	0	2.589	−1.120 × 10^2	−3.012 × 10^−1

---

Fig. 8 compares the fitted monodisperse PhaseTime curves with experimental phase boundaries for PVME-based blends. Fig. 8a shows spinodal data for two poly(vinyl methyl ether)/deuterated polystyrene blends, denoted PVME_PSD_L and PVME_PSD_H, which differ in molecular weight and polydispersity;^38,68 the L system is characterized by lower molecular weights and narrower molecular-weight distributions, whereas the H system involves higher molecular weights and broader distributions, particularly in PSD.³⁸ The effective monodisperse representation captures the location and curvature of both reported spinodal boundaries. Fig. 8b shows the corresponding fits to binodal data for poly(vinyl methyl ether)/polystyrene blends with two different polystyrene molecular weights, PVME_PS62k and PVME_PS100k.⁶⁹ In both cases, the fitted curves closely follow the experimental coexistence boundaries.

Fig. 8 Experimental phase boundaries and corresponding monodisperse PhaseTime fits for PVME-based blends: (a) spinodal curves of PVME_PSD_L and PVME_PSD_H; (b) binodal curves of PVME_PS62k and PVME_PS100k. Fitted model parameters are listed in Table 3.

Fig. 9 shows the fitted PhaseTime spinodal curve for the deuterated polybutadiene/protonated polybutadiene blend PBD_DPBD, which exhibits UCST behaviour.⁷⁰ The fitted monodisperse model captures the reported spinodal boundary for this chemically distinct blend.

Fig. 9 Experimental spinodal data and corresponding PhaseTime fit for the deuterated polybutadiene/protonated polybutadiene blend PBD_DPBD, exhibiting UCST behaviour. Fitted model parameters are listed in Table 3.

Fig. 10 shows the phase diagrams of acetone/polystyrene blends with two different polystyrene molecular weights.⁷¹ This system displays more complex phase behaviour that depends strongly on the chain size of polystyrene. For the lower-molecular-weight case, acetone_PS5k, the phase diagram exhibits combined LCST–UCST behaviour; see Fig. 10a. Increasing the polystyrene molecular weight to approximately 20 kDa changes the topology to an hourglass-shaped phase diagram. The fitted PhaseTime curves reproduce both the change in topology and the main features of the experimental binodal boundaries.

Fig. 10 Experimental binodal curves and corresponding PhaseTime fits for acetone/polystyrene blends with different polystyrene molecular weights: (a) acetone_PS5k, showing combined LCST-UCST behaviour; (b) acetone_PS20k, showing hourglass behaviour. Fitted model parameters are listed in Table 3.

Fig. 11 shows the fitted PhaseTime binodal curve for the water/PVA system, which exhibits closed-loop phase behaviour.⁷² Compared with the previous examples, this system presents a more intricate coexistence topology, with both upper and lower bounds of immiscibility. The fitted model captures the closed-loop boundary, demonstrating that the same thermodynamic and numerical framework can represent nonmonotonic phase behaviour across chemically diverse systems.

Fig. 11 Experimental binodal curve and corresponding PhaseTime fit for the water/PVA system, exhibiting closed-loop phase behaviour. Fitted model parameters are listed in Table 3.

Overall, these results show that the effective monodisperse formulation can reproduce a wide range of experimental phase-boundary topologies with high accuracy after case-specific fitting. The fitted coefficients, however, should be viewed as system-dependent effective parameters rather than transferable descriptors, since the effects of molecular-weight distribution are absorbed implicitly into the fit. The monodisperse description thus provides only a reduced representation of the molecular constitution, whereas an explicit polydisperse treatment is expected to be more transferable across systems with different dispersities and is essential when fractionation, cloud-shadow asymmetry, or reaction-induced evolution of the molecular-weight distribution must be described.⁷

3.2 Phase behavior under reaction conditions

3.2.1 Monodisperse approximation. We first examine reaction-induced phase behaviour within an effective monodisperse description, in which the evolving molecular distributions are represented by reaction-dependent effective chain sizes. For validation, we consider the epoxy/rubber system discussed in the review by Williams et al.,⁴ originally based on the work cited therein.³

The system consists of a DGEBA-based diepoxide, a stoichiometric amount of the diamine hardener (3DCM), and a rubber phase based on a statistical copolymer of butadiene and acrylonitrile. In the present representation, the rubber is assigned a fixed chain size of N₂ = 13.73, while the interaction parameter is taken as χ = 0.63 ± 0.03 at 348 K, as reported from fitting cloud-point conversion data.³ The evolution of the epoxy chain size N₁(a₁) is described by the Carothers-type relation in eqn (23), using γ_c = 4/3.

The corresponding phase diagram is shown in Fig. 12 in composition–reaction coordinates. At low extent of reaction a₁, the blend remains in the single-phase region, whereas phase separation is induced as the reaction proceeds and the effective molecular weight of the epoxy increases. According to the present calculation, the critical point is located at (φ₂, a₁) = (0.24, 0.205), in agreement with the value reported by Williams et al.⁴

Fig. 12 Calculated phase diagram for the DGEBA/3DCM/rubber system under reaction conditions, shown in composition–reaction coordinates using the effective monodisperse approximation. For details see text.

3.2.2 Polydisperse systems with Flory–Stockmayer distributions. Having examined reaction-induced phase behaviour within the effective monodisperse approximation, we now consider explicitly polydisperse blends, in which the evolving molecular constitution of the reactive species is described by the Flory–Stockmayer model introduced in section 2.4.2. The aim is to assess how the broadening of the molecular-weight distribution under reaction affects the predicted phase boundaries, and to evaluate the role of distribution truncation and moment-preserving corrections in the numerical treatment.

As a representative test case, species 2 is taken to be monodisperse with degree of polymerization N₂ = 20, whereas species 1 evolves according to a Flory–Stockmayer chain-size distribution. At the initial state, a₁(t = 0) = 0, species 1 consists of two monomeric precursor classes: A-type molecules with relative concentration A_fi = 0.4, functionality f_i = 3, and reduced molecular volume N^A_i = 0.5; and B-type molecules with relative concentration B_gj = 0.6, functionality g_j = 2, and reduced molecular volume N^B_j = 0.8, whereas constituents with relative volume-fraction contributions below w_min = 10⁻⁹ were omitted; for a detailed description refer to section S3.3 (SI). Throughout the protocol, the temperature is held constant at T = 450 K, while the extent of reaction increases linearly according to a₁(t) = 0.4t/100. Here, χ is composition-independent, and follows eqn (21) with d₀ = 0.76, d₁ = −117 K and d₂ = 0 yielding LCST behavior.

As the reaction proceeds, the Flory–Stockmayer distribution broadens substantially and, at the gel point, develops an infinite tail. In practice, the discrete distribution must therefore be truncated at a finite maximum chain size. If the cutoff is too low, omission of the long-chain tail may distort the predicted phase behaviour. To reduce this error, a moment-matching correction is introduced through additional effective chain lengths, such that the moments N_1n, N_1w, and N_1z of the truncated distribution are restored to their analytical Flory–Stockmayer values. Details of this procedure are given in sections S3.6–S3.9 (SI).

The resulting phase diagram is shown in Fig. 13 in reaction-composition coordinates. Panel (a) shows the cloud curves together with the corresponding spinodal, critical line, and critical point, while panel (b) shows the associated shadow curves. Results are presented for multiplicity values of n_max = m_max = 5, 10, 20, and 40, with moment matching applied in all cases.

Fig. 13 Reaction-composition phase diagram for the representative polydisperse reactive blend described by a Flory–Stockmayer distribution for species 1 and a monodisperse species 2 with N₂ = 20. (a) Cloud curves for different cutoff values of the maximum chain size, n_max = m_max = 5, 10, 20, and 40, shown together with the corresponding spinodal, critical line, critical point, and the binodal of the monodisperse counterpart. (b) Corresponding shadow curves for the same cutoff values. In all cases, the moment-matching correction was applied so that the moments N_1n, N_1w, and N_1z of the truncated distributions coincide with their analytical values.

As the cutoff increases, the cloud and shadow curves broaden and approach convergence, whereas smaller cutoffs lead to progressively narrower coexistence regions. The effect is most evident on the right branch of the cloud curve, while the left branch is only minimally affected. By contrast, the shadow curves change only weakly with cutoff. In particular, the left shadow branch, although expected to be more sensitive, occurs at very low volume fractions, so the corresponding differences are barely visible. For comparison, the corresponding monodisperse binodal is significantly narrower, highlighting the effect of explicit polydispersity on the predicted phase boundary.

Because the moment-matching correction restores the N_1n, N_1w, and N_1z exactly, the spinodal and critical quantities are unchanged by the cutoff and therefore collapse onto the same curves in Fig. 13. This is consistent with the analysis of sections 2.5.1 and 2.5.2, where these quantities depend only on the corresponding distribution moments.

In contrast, the cloud and shadow curves remain sensitive to the detailed form of the truncated distribution, particularly at low cutoff values. Illustrative results without moment matching are given in section S5 (SI); in that case, the spinodal, critical point, and corresponding effective monodisperse binodal also become cutoff-dependent, and the approximation deteriorates much more strongly as the cutoff is reduced. The shadow curves are also affected appreciably as well, unlike in the moment-matched case where their cutoff dependence is weak.

4. Concluding remarks

This work introduced PhaseTime, a Python framework for computing time- or protocol-dependent phase diagrams of reactive polydisperse polymer blends within a Flory–Huggins-type thermodynamic description. The framework combines quasibinary thermodynamics, composition- and temperature-dependent interaction models, evolving molecular-distribution functions, and numerical procedures for evaluating spinodal curves, critical points, binodal curves, and cloud/shadow curves. The protocol coordinate may represent time, temperature, conversion, or another externally prescribed variable, allowing molecular constitution and thermodynamic parameters to be updated consistently along a reaction or processing path.

The examples demonstrate that the framework can reproduce several idealized phase-diagram topologies, including UCST, LCST, combined UCST/LCST, closed-loop, and hourglass-type behavior. The benchmark fitting cases show how interaction-model parameters can be estimated from literature phase-boundary data. Under reaction conditions, the calculations illustrate how evolving chain size, molecular-weight distribution, and interaction parameters shift stability and coexistence boundaries. Comparisons between monodisperse approximations and polydisperse Flory–Stockmayer distributions further show that molecular-distribution effects can influence the predicted phase behavior, particularly in cloud/shadow calculations.

Beyond the thermodynamic formulation, PhaseTime provides a reproducible computational workflow through a common backend exposed by both a command-line interface and a Python API. The implementation includes continuation procedures, curve segmentation, diagnostic grid calculations, plotting utilities, and parameter-estimation workflows. These tools support the analysis of phase diagrams with multiple branches, disconnected segments, hidden coexistence curves, or complex critical-point structures, and make the framework useful for comparing model assumptions and inspecting numerical solutions.

The present implementation remains within a mean-field Flory–Huggins framework and therefore neglects fluctuation corrections and renormalization effects, which may become important near critical points or in strongly asymmetric disperse blends. Such effects have been treated separately using renormalized descriptions of polymer blends, but are outside the scope of the present work.^79,80 Moreover, the computed phase diagrams should be interpreted as quasi-equilibrium thermodynamic references, corresponding to the limiting case in which demixing and diffusive relaxation are fast compared with the imposed protocol. They do not describe spatial morphology evolution, domain coarsening, nucleation barriers, vitrification, or gelation kinetics directly.

Future work will first focus on developing a phase-field model for quasibinary blends under steady-state conditions, using PhaseTime phase diagrams as reference solutions for the thermodynamic limit. This model will then be extended to reaction-induced phase separation by coupling the spatial evolution to reaction kinetics, evolving molecular distributions, and protocol-dependent interaction parameters. In this setting, the framework will provide the quasi-equilibrium reference for benchmarking and interpreting the spatial simulations.

Author contributions

Aristotelis P. Sgouros: conceptualization, methodology, software, validation, resources, formal analysis, investigation, data curation, writing – original draft, writing – review & editing, visualization, supervision, project administration, funding acquisition. Konstantinos S. Gkoutis: software, validation, writing – review & editing. Anthony Bocahut: validation, writing – review & editing, project administration. Eléonore Mathis: validation, writing – review & editing. Doros N. Theodorou: conceptualization, methodology, validation, resources, writing – original draft, writing – review & editing, supervision, project administration, funding acquisition.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

Supplementary information is available. See DOI: https://doi.org/10.1039/d6me00077k.

The supplementary information (SI) contains a pdf file with the detailed theoretical and numerical material supporting this work. Section S1 presents the derivation of the chemical potentials for the polydisperse constituents 1i and 2j. Section S2 derives the cloud-point equations used to compute the cloud and shadow curves. Section S3 describes the Flory–Stockmayer distribution implemented in PhaseTime. Section S4 reports the curve-segmentation algorithm used for plotting and parameter fitting. Section S5 provides additional reaction–composition phase diagrams computed without moment matching. Ref. 6, 7, 33, 38, 40, 41, 49–51 and 81 are cited in the SI. In addition, the code, input files, and data required to reproduce the calculations reported in this article are provided in the SI as a zipped frozen snapshot of the PhaseTime GitHub repository;⁴⁷ the snapshot corresponds to Git tag v0.2-submission-msde and commit 05899f75264715e782371cbd62f77a953a319418. Installation instructions, tests, and example workflows are provided in the repository README files.

Acknowledgements

A. P. Sgouros thanks Solvay Specialty Polymers USA, for financial support. A. P. Sgouros is grateful to Dr. A. Bocahut and Dr. E. Mathis during his productive stay at SYENSQO in March 2025. A. P. Sgouros gratefully acknowledges Dr. P. Dallas and Dr. E. Sakellis for useful discussions and input on the experimental aspects of this work. PhaseTime was tested on the Helios HPC cluster in the framework of the project “Center of Excellence for Advanced Organic Materials, Theranostics in Cancer and Bioelectronics” (Project ID: 5173962). AI-based tools were used to assist with code refinement and language editing. All algorithms and results were independently developed, verified, and validated by the authors.

References

1. A. R. Ajitha and S. Thomas, in Compatibilization of Polymer Blends, ed. A. R. Ajitha and S. Thomas, Elsevier Inc., 2020, pp. 1–29.
2. Z. Sun, L. Xu, Z. Chen, Y. Wang, R. Tusiime, C. Cheng, S. Zhou, Y. Liu, M. Yu and H. Zhang, Polymers, 2019, 11, 461.
3. S. M. Moschiar, C. C. Riccardi, R. J. J. J. Williams, D. Verchere, H. Sautereau and J. P. Pascault, J. Appl. Polym. Sci., 1991, 42, 717–735.
4. R. J. J. Williams, B. A. Rozenberg and J. P. Pascault, Adv. Polym. Sci., 1997, 128, 95–156.
5. S. J. Mumby, P. Sher and B. E. Eichinger, Polymer, 1993, 34, 2540–2545.
6. S. J. Mumby and P. Sher, Macromolecules, 1994, 27, 689–694.
7. S. J. Mumby, P. Sher and J. van Ruiten, Polymer, 1995, 36, 2921–2927.
8. N. Clarke, T. C. B. McLeish and S. D. Jenkins, Macromolecules, 1995, 28, 4650–4659.
9. Q. Guo and Z. Liu, J. Therm. Anal. Calorim., 2000, 59, 101–120.
10. B. J. P. Jansen, K. Y. Tamminga, H. E. H. Meijer and P. J. Lemstra, Polymer, 1999, 40, 5601–5607.
11. S. L. Veatch and S. L. Keller, Biochim. Biophys. Acta, Mol. Cell Res., 2005, 1746, 172–185.
12. M. Seiler, Fluid Phase Equilib., 2006, 241, 155–174.
13. N. Z. Tomić, A. D. Marinković, Đ. Veljović, K. Trifković, S. Lević, V. Radojević and R. J. Heinemann, Int. J. Adhes. Adhes., 2018, 81, 11–20.
14. D. J. Duffy, A. M. Heintz, H. D. Stidham, S. L. Hsu, W. Suen, W. Chu and C. W. Paul, J. Adhes., 2003, 79, 1091–1107.
15. F. Siepmann, J. Siepmann, M. Walther, R. J. MacRae and R. Bodmeier, J. Controlled Release, 2008, 125, 1–15.
16. J. B. Nephew, T. C. Nihei and S. A. Carter, Phys. Rev. Lett., 1998, 80, 3276–3279.
17. E. P. Favvas and A. C. Mitropoulos, J. Eng. Sci. Technol. Rev., 2008, 1, 25–27.
18. K. Binder, J. Chem. Phys., 1983, 79, 6387–6409.
19. J. S. Langer, Ann. Phys., 1971, 65, 53–86.
20. I. M. Lifshitz and V. V. Slyozov, J. Phys. Chem. Solids, 1961, 19, 35–50.
21. C. Wagner, Zeitschrift für Elektrochemie, Berichte der Bunsengesellschaft für physikalische Chemie, 1961, 65, 581–591.
22. W. Li, A. J. Ryan and I. K. Meier, Macromolecules, 2002, 35, 5034–5042.
23. R. W. Venderbosch, H. E. H. Meijer and P. J. Lemstra, Polymer, 1994, 35, 4349–4357.
24. S. Swier and B. Van Mele, Polymer, 2003, 44, 2689–2699.
25. E. Girard-Reydet, H. Sautereau, J. P. Pascault, P. Keates, P. Navard, G. Thollet and G. Vigier, Polymer, 1998, 39, 2269–2279.
26. P. J. Flory, J. Chem. Phys., 1942, 10, 51–61.
27. M. L. Huggins, J. Chem. Phys., 1941, 9, 440.
28. D. Qian, T. C. T. Michaels and T. P. J. Knowles, J. Phys. Chem. Lett., 2022, 13, 7853–7860.
29. J. P. de Souza and H. A. Stone, J. Chem. Phys., 2024, 161, 44902.
30. K. Kamide and Y. Miyazaki, Polym. J., 1981, 13, 325–341.
31. K. Kamide, T. Abe and Y. Miyazaki, Polym. J., 1982, 14, 355–361.
32. K. Kamide, S. Matsuda, T. Dobashi and M. Kaneko, Polym. J., 1984, 16, 839–855.
33. R. Koningsveld and A. J. Staverman, J. Polym. Sci., Polym. Phys. Ed., 1968, 6, 325–347.
34. K. Šolc, Macromolecules, 1970, 3, 665–673.
35. 3PDB Flory–Huggins, https://pppdb.uchicago.edu/flory, (accessed 29 April 2026).
36. Flory Huggins Phase Diagram Of Polymer Mixing, https://demonstrations.wolfram.com/FloryHugginsPhaseDiagramOfPolymerMixing/, (accessed 29 April 2026).
37. Y. Qiang and D. Zwicker, J. Open Source Softw., 2025, 10, 7388.
38. C. Qian, S. J. Mumby and B. E. Eichinger, Macromolecules, 1991, 24, 1655–1661.
39. W. H. Carothers, Trans. Faraday Soc., 1936, 32, 39–49.
40. W. H. Stockmayer, J. Polym. Sci., 1952, 9, 69–71.
41. W. H. Stockmayer, J. Polym. Sci., 1953, 11, 424.
42. P. K. Chan and A. D. Rey, Macromolecules, 1997, 30, 2135–2143.
43. T. Higuchi, Y. Yano, T. Aita, S. Takami and T. Adschiri, J. Chem. Eng. Jpn., 2013, 46, 709–715.
44. P. K. Chan and A. D. Rey, Macromolecules, 1996, 29, 8934–8941.
45. J. C. Lee, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 1930–1935.
46. C. Li and A. Strachan, Polymer, 2018, 149, 30–38.
47. A. P. Sgouros, PhaseTime, https://github.com/ArisSgouros/PhaseTime, (accessed 21 June 2026).
48. P. J. Flory, Principles of polymer chemistry, Cornell University Press, 1953.
49. C. W. Macosko and D. R. Miller, Macromolecules, 1976, 9, 199–206.
50. R. Odle, J. D. Mitchell and J. T. Bendler, J. Polym. Sci., Part B:Polym. Phys., 2019, 57, 1415–1422.
51. R. Bachmann and J. T. Bendler, Macromol. Theory Simul., 2025, 34, 2400073.
52. J. Sanders, Veusz, version 3.6.2, scientific plotting software. Available at: https://veusz.github.io/, (accessed 20 June, 2026).
53. I. Kutskir, Photopea, https://www.photopea.com/, (accessed 29 April 2026).
54. R. P. Brent, Algorithms for minimization without derivatives, Courier Corporation, 2013.
55. P. Virtanen, R. Gommers, T. E. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, A. Vijaykumar, A. Pietro Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Häggström, C. Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G.-L. Ingold, G. E. Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavič, J. Nothman, J. Buchner, J. Kulick, J. L. Schönberger, J. V. de Miranda Cardoso, J. Reimer, J. Harrington, J. L. C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. Kümmerer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y. O. Halchenko, Y. Vázquez-Baeza and SciPy 1.0 Contributors, Nat. Methods, 2020, 17, 261–272.
56. M. A. Branch, T. F. Coleman and Y. Li, SIAM J. Sci. Comput., 1999, 21, 1–23.
57. D. J. Kozuch, W. Zhang and S. T. Milner, Polymers, 2016, 8, 241.
58. A. Milchev, S. A. Egorov, J. Midya, K. Binder and A. Nikoubashman, ACS Macro Lett., 2020, 9, 1779–1784.
59. E. Patyukova, E. Xi and M. R. Wilson, Macromolecules, 2021, 54, 2763–2773.
60. H. Jacobson, C. O. Beckmann and W. H. Stockmayer, J. Chem. Phys., 1950, 18, 1607–1612.
61. M. Lang and K. S. Kumar, Macromolecules, 2021, 54, 7021–7035.
62. P. G. De Gennes, J. Chem. Phys., 1980, 72, 4756–4763.
63. H. H. Kausch and M. Tirrell, Annu. Rev. Mater. Sci., 1989, 19, 341–377.
64. R. S. Qin and H. K. Bhadeshia, Mater. Sci. Technol., 2010, 26, 803–811.
65. S.-P. Lyu, J. J. Cernohous, F. S. Bates and C. W. Macosko, Macromolecules, 1999, 32, 106–110.
66. N. Hansen and A. Ostermeier, Evol. Comput., 2001, 9, 159–195.
67. N. Hansen, The CMA Evolution Strategy: A Tutorial, arXiv, 2023, preprint, arXiv:1604.00772v2,  DOI:10.48550/arXiv.1604.00772.
68. C. C. Han, B. J. Bauer, J. C. Clark, Y. Muroga, Y. Matsushita, M. Okada, Q. Tran-Cong, T. Chang and I. C. Sanchez, Polymer, 1988, 29, 2002–2014.
69. D. J. Walsh, G. T. Dee, J. L. Halary, J. M. Ubiche, M. Millequant, J. Lesec and L. Monnerie, Macromolecules, 1989, 22, 3395–3399.
70. S. Sakurai, H. Hasegawa, T. Hashimoto, I. G. Hargis, S. L. Aggarwal and C. C. Han, Macromolecules, 1990, 23, 451–459.
71. K. S. Siow, G. Delmas and D. Patterson, Macromolecules, 1972, 5, 29–34.
72. H.-G. Elias, Macromolecules. Volume 1: Structure and Properties, Plenum Press, New York, 1984.
73. T. Shiomi, F. Hamada, T. Nasako, K. Yoneda, K. Imai and A. Nakajima, Macromolecules, 1990, 23, 229–233.
74. H. Höcker, G. J. Blake and P. J. Flory, Trans. Faraday Soc., 1971, 67, 2251–2257.
75. Agency for Toxic Substances and Disease Registry (US), Physical and Chemical Properties of Acetone, https://www.ncbi.nlm.nih.gov/books/NBK590394/table/ch4.tab2/, (accessed 29 April 2026).
76. National Institute of Standards and Technology, Composition of polystyrene, https://pml.nist.gov/cgi-bin/Star/compos.pl?matno=226&mode=text&refer=ap, (accessed 29 April 2026).
77. National Institute of Standards and Technology, Composition of polyvinyl alcohol, https://pml.nist.gov/cgi-bin/Star/compos.pl?matno=230, (accessed 29 April 2026).
78. Softonic, CurveSnap for Windows, https://curvesnap.en.softonic.com/, (accessed 29 April 2026).
79. P. M. Rauscher, J. Chem. Phys., 2023, 159, 244906.
80. J. Qin and D. C. Morse, J. Chem. Phys., 2009, 130, 224902.
81. B. H. Zimm and W. H. Stockmayer, J. Chem. Phys., 1949, 17, 1301–1314.

---