Simulation of Binodal and Spinodal Curves of Phase State Diagrams for Binary Polymer Systems

A new approach is proposed for simulating binodal and spinodal curves of phase diagrams for binary polymer systems. It is shown that the Flory–Huggins theory makes it possible to predict phase behavior in a wide range of temperatures and concentrations based on limited data on the components’ solubility. The approbation data of the technique are presented in the example of PS–PB and PS–PMMA systems, for which generalized phase diagrams are constructed.


Introduction
The construction of phase state diagrams has been and remains a priority task both in the thermodynamics of polymer solutions and melts and in polymer materials science [1][2][3][4][5]. It is known that the boundary curves, liquidus, and solidus lines are of fundamental importance for the construction and interpretation of phase state diagrams since they determine not only the type of phase equilibrium but also the arias of formation of various phase structures in the temperature-concentration regions.
Thus, the binodal curve in systems with amorphous equilibrium separates the region of the true solutions from the region of the two-phase state. The region located between the binodal curve and the spinodal curve makes it possible to identify the position on the diagram of the region of the metastable and labile solutions. Finally, in the middle region of the diagram, there are boundary curves framing the region of phase reversal. The construction of generalized diagrams [6] makes it possible to determine the regions of heterogeneous fluctuations, the position of critical points, and the regions of glass transition, viscous flow, and the thermal stability of components. Thus, with a generalized phase state diagram, it is possible to predict the structural and morphological transformations that occur under certain conditions of processing (high-temperature regimes with rapid cooling), the exploitation of materials, and also during the synthesis of copolymers. It is obvious that to solve some material science problems, it is necessary to expand the temperature and concentration intervals for plotting state diagrams.
The importance of predicting and constructing state diagrams of polymer-polymer systems is confirmed by the fact that two independent directions are currently being developed in this area of the physical chemistry of polymers. The first-theoretical simulation-concerns the statistical thermodynamics of solutions and polymer blends. These are the classical Flory-Huggins theory, the new Flory theory, the Sanchez-Lacombe theory, and Prigogine-Paterson theory [1][2][3][4][5], the statistical associating fluid theory (SAFT) [7], the polymer reference interaction site model (PRISM) [8], and the lattice cluster theory (LCT) [9]. Recently, the adapted mean field theory [10], the hybrid molecular dynamics-Monte Carlo method [11], the self-consistent field theory (SCFT) [12], and theoretically informed Langevin dynamics [13] have been used in diagram simulation. All of these theories make it possible to obtain analytical dependences of the free energy of mixing components on the composition of systems and to determine the position of binodal and spinodal boundary curves and critical states in the framework of classical thermodynamics. It is interesting to note that comparative studies of the prediction of the above statistical theories of binary polymer compatibility have shown that the classical Flory-Huggins theory has a good ability to predict phase equilibrium in polymer solutions and melts [2].
The second direction is associated with the independent determination of the Flory-Huggins parameter based on the data of inverse gas chromatography, polymer solutions, intrinsic viscosity, and sorption [3,[14][15][16][17]. That is, the prediction of the boundary curves here is partly based on preliminary experimental data for determining the interaction parameter. Extrapolation of the obtained physicochemical characteristics to the temperature ranges of interest to researchers also implies the calculation of the free energy of mixing and the subsequent expansion of the investigated range of the components' compatibility. It is obvious that the construction of such diagrams is identical to the prediction of segments of phase-state diagrams. In this case, the biggest difficulties arise with the construction of binodal curves due to the uncertainty in estimating the compositions of coexisting phases using binodal equations.
In this work, we chose the second approach, but we estimated the interaction parameters not by third-party methods, but from the data on the solubility of the components in a limited range. At the same time, we propose a more efficient and convenient method for numerically determining the compositions of coexisting phases and constructing the position of binodal and spinodal curves in a wide range of temperatures and molecular weights of binary system components. This approach makes it possible to construct generalized phase diagrams for different types of polymer-polymer systems, which has been tested on a number of examples.

Theoretical Methodology
The Flory-Huggins theory of polymer solutions [1][2][3][4][5] is based on the simplest lattice model of an athermal solution, which contains assumptions about the volumes of macromolecular chain segments and their relative mobility, which is associated with the combinatorial entropy of all possible permutations within such a lattice.
Traditionally, in the framework of the Flory-Huggins theory, the free energy of mixing two polymers is represented as follows: where R is the universal gas constant, T is the temperature, V is the polymer mixing volume, V r is the reference volume, usually taken as the molar volume of a repeating unit in the system (in calculations, a reference volume of 100 cm 3 /mol is traditionally used [2]), ϕ 1 and ϕ 2 are the volume fractions of the first and second components in the system, r 1 and r 2 are their polymerization degrees, χ 12 is the Flory-Huggins interaction parameter.
Expressions for chemical potentials are as follows: The conditions at the critical point are as follows: Spinodal concentrations are defined as the points at which the second derivative of the free energy of mixing with respect to concentration is zero. The spinodal equation takes the form: 1 where χ is the averaged value of the pair interaction parameter. The binodal equations are calculated from the condition that the chemical potentials of both polymers are equal in two coexisting phases: where and " refer to different coexisting phases. For thermodynamic analysis of the obtained experimental data on binodal curves (or their fragments), the presented equations lead to the following expressions for the pair interaction parameter: The averaged value of the pair interaction parameter, assuming the absence of its concentration dependence, can be represented as follows: Figure 1 shows a scheme of the methodology for the thermodynamic analysis of binary polymer systems and the construction of generalized phase diagrams. At the initial stage, fragments of binodal curves are plotted based on the experimentally obtained data on the solubility of polymers in each other ( Figure 1a). Then, according to the compositions of the coexisting phases ϕ 1 and ϕ 1 at different temperatures T, using Equation (9), the numerical values of the pair interaction parameter χ are calculated, and its temperature dependence is plotted (Figure 1b). Extrapolation of this dependence to χ cr makes it possible to obtain information about the critical solution temperature of the components (in this case, the upper critical solution temperature -UCST). In coordinates χ = A + B· 1 T it is possible to identify χ data in a wide range of temperatures at T ≥ T ex , T ≤ T ex , T ≥ USCT. Based on these data, using Equations (5) and (6), the boundary spinodal (dashed line) and binodal (solid line) curves are calculated, and a generalized phase diagram is constructed (Figure 1c) in the selected temperature range, which traditionally denotes the region of homogeneous states (I), the region of labile structures (II), and the region of metastable states (III).
We should consider in more detail the procedure for determining data on the coexisting phase compositions. It is known that the concentration dependence of free energy for a polymer-polymer system at a temperature T ≤ UCST, corresponding to the limited mixing of components, is a complex curve with two minima and two inflection points ( Figure 2a). By drawing a common tangent to the dependence ∆G − ϕ, corresponding to the equality of chemical potentials, it can determine the concentrations at the points of the curves contact (points 2) about thus assess the information about the coexisting phases compositions for the binodal curve (ϕ b and ϕ b ), and in inflection points (points 3) receive information for constructing a spinodal curve (ϕ c and ϕ c ). Numerical determination of binodal and spinodal concentrations of more or less symmetric concentration dependences of the free mixing energy does not cause any difficulties [18]. This situation is realized in cases where the molecular weights of the components do not differ much. When analyzing a similar dependence for a binary system with significantly different molecular weights of components-for example, a polymer solution, a polymer-plasticizer mixture, a polymer-oligomer system, and even a high-molecular polymer-low-molecular polymer blend-significant difficulties arise in the standard approach. It can be seen (Figure 2b) that in this case the dependence ∆G m (ϕ) becomes sharply asymmetric, one of the binodal concentrations shifts to the region of small values, the curve flattens out in this place, which greatly complicates the numerical determination of the phase composition. In this case, the direct mathematical solution of the system of Equation (6) leads to a multiplicity of solutions that are difficult to automatically limit. We should consider in more detail the procedure for determining data on the coexisting phase compositions. It is known that the concentration dependence of free energy for a polymer-polymer system at a temperature T ≤ UCST, corresponding to the limited mixing of components, is a complex curve with two minima and two inflection points ( Figure 2a). By drawing a common tangent to the dependence Δ − , corresponding to the equality of chemical potentials, it can determine the concentrations at the points of the curves contact (points 2) about thus assess the information about the coexisting phases compositions for the binodal curve ( and ), and in inflection points (points 3) receive information for constructing a spinodal curve ( and ). Numerical determination of binodal and spinodal concentrations of more or less symmetric concentration dependences of the free mixing energy does not cause any difficulties [18]. This situation is realized in cases where the molecular weights of the components do not differ much. When analyzing a similar dependence for a binary system with significantly different molecular weights of components-for example, a polymer solution, a polymer-plasticizer mixture, a polymer-oligomer system, and even a high-molecular polymer-low-molecular polymer blend-significant difficulties arise in the standard approach. It can be seen (Figure 2b) that in this case the dependence Δ becomes sharply asymmetric, one of the binodal concentrations shifts to the region of small values, the curve flattens out in this place, which greatly complicates the numerical determination of the phase composition. In this case, the direct mathematical solution of the system of Equation (6) leads to a multiplicity of solutions that are difficult to automatically limit.  We should consider in more detail the procedure for determining data on the coexisting phase compositions. It is known that the concentration dependence of free energy for a polymer-polymer system at a temperature T ≤ UCST, corresponding to the limited mixing of components, is a complex curve with two minima and two inflection points (Figure 2a). By drawing a common tangent to the dependence Δ − , corresponding to the equality of chemical potentials, it can determine the concentrations at the points of the curves contact (points 2) about thus assess the information about the coexisting phases compositions for the binodal curve ( and ), and in inflection points (points 3) receive information for constructing a spinodal curve ( and ). Numerical determination of binodal and spinodal concentrations of more or less symmetric concentration dependences of the free mixing energy does not cause any difficulties [18]. This situation is realized in cases where the molecular weights of the components do not differ much. When analyzing a similar dependence for a binary system with significantly different molecular weights of components-for example, a polymer solution, a polymer-plasticizer mixture, a polymer-oligomer system, and even a high-molecular polymer-low-molecular polymer blend-significant difficulties arise in the standard approach. It can be seen (Figure 2b) that in this case the dependence Δ becomes sharply asymmetric, one of the binodal concentrations shifts to the region of small values, the curve flattens out in this place, which greatly complicates the numerical determination of the phase composition. In this case, the direct mathematical solution of the system of Equation (6) leads to a multiplicity of solutions that are difficult to automatically limit. Analyzing the principle of equality of chemical potentials in different phases and its interpretation as a common tangent to the dependence ∆G m (ϕ), we can assume that there are the same values at two points of the concentration dependences of chemical potentials. This means that the dependence ∆µ 2 (∆µ 1 ) must have a self-intersection point (two points of the dependence must have the same value). For a single-phase system, this dependence has no self-intersection points, and for a two-phase system, it has a single point corresponding to the equality of chemical potentials in different phases, i.e., thermodynamic condition of phase equilibrium. Figure 3 shows such a schematic concentration dependence interpretation as a common tangent to the dependence Δ , we can assume that there are the same values at two points of the concentration dependences of chemical potentials. This means that the dependence Δ Δ must have a self-intersection point (two points of the dependence must have the same value). For a single-phase system, this dependence has no self-intersection points, and for a two-phase system, it has a single point corresponding to the equality of chemical potentials in different phases, i.e., thermodynamic condition of phase equilibrium. Figure 3 shows such a schematic concentration dependence of chemical potentials where the values of Δ and Δ correspond to the same value of the concentration of the components. Therefore, the enumeration of concentrations in the range from 0 to 1 gives a set of pairs of chemical potentials for this dependence. The point of self-intersection (1) determines, in this case, the concentrations of the binodal curve. Accordingly, when constructing a generalized phase diagram, the spinodal Equation (5) can be used to determine the composition of the coexisting phases of the spinodal curve. The transformation of this equation allows us to express its roots in the following form: To find the binodal points, we propose using the following mathematical algorithm:  Accordingly, when constructing a generalized phase diagram, the spinodal Equation (5) can be used to determine the composition of the coexisting phases of the spinodal curve. The transformation of this equation allows us to express its roots in the following form: To find the binodal points, we propose using the following mathematical algorithm:  Figure 1b).

Experiment
The following monodisperse polymers were used as objects for approbation of the methodology: polystyrene (PS) with an average molecular weight of M w = 37.8 kDa and M w = 0.8 kDa (PS-1 and PS-2, respectively) manufactured by Waters Associates, polybutadiene (PB) with M w = 61.6 kDa (Aldrich), and polymethyl methacrylate (PMMA) with M w = 89 kDa (Glass).
The solubility data of the components were obtained using the optical interferometry method [19]. Pressed polymer films about 5 × 10 mm in size were placed between two optically transparent glasses, on the inside of which a translucent metal layer (nichrome) was deposited by thermal vacuum deposition. Glasses with polymers between them were clamped in a temperature-controlled cell (±1 K) at a temperature slightly above the glass transition (melting) temperature of polymers to a variable thickness of 100-120 µm, fixed by special metal wedges. The polymers were in optical contact with the inner surfaces of the glasses. A digital camera through a microscope with standard magnifications fixed the interference from a light source (helium-neon laser with a wavelength of 632 nm), passing perpendicular to the film contact plane and forming a concentration profile of polymer interaction at the contact boundary. The measurements were carried out in the mode of a stepwise increase and decrease in temperature from 400 to 540 K. The temperature step was 20-40 K with holding at each step for at least 30-60 min, which makes it possible to obtain equilibrium data. The reproducibility of data in heating-cooling cycles allows us to speak about the equilibrium values of the compositions of coexisting phases at each temperature. This technique makes it possible to simultaneously observe the entire concentration range at the contact boundary of the two components, and the error in determining the compositions of the coexisting phases at each temperature is about 5%. The procedure for conducting the experiment and processing the interferograms did not differ from the traditional procedure [20][21][22][23].

Results and Discussion
To test the calculations using the proposed method, experimental results on solubility in the PS-PB and PS-PMMA systems were used. Figure 4 shows fragments of binodal curves for the PB-PS-1 and PMMA-PS-2 systems. The molecular weights of the first pair do not differ much (M w PB /M w PS-1 = 1.6), but in the second case, they differ significantly (M w PMMA /M w PS-2 = 111.3). Thus, we tried to simulate two different cases for polymer-polymer systems according to the scheme in Figure 2. The data for both systems are characterized by limited solubility, in the range from 400 to 530 K for the PB-PS-1 system and from 460 to 520 K for the PMMA-PS-2 system. For both pairs, the experimental range of values was limited by the glass transition/melting temperatures from below and the thermal oxidative degradation from above. Nevertheless, it can be said that the solubility of the components increases with temperature growth, and it should expect the appearance of an upper critical solution temperature (UCST) in the higher temperatures' region.  According to the method proposed above, using Equation (9), pair interaction parameters were calculated, and their temperature dependences were plotted ( Figure 5). It can be seen that both systems demonstrate a linear change χ in the proposed coordinates, and interpolation of such a dependence is possible in the form = 3.4748 • According to the method proposed above, using Equation (9), pair interaction parameters were calculated, and their temperature dependences were plotted ( Figure 5). It can be seen that both systems demonstrate a linear change χ in the proposed coordinates, and interpolation of such a dependence is possible in the form χ = 3.4748·(1/T) − 0.0022 for the PB-PS-1 system and χ = 142.78·(1/T) − 0.1608 for the PMMA-PS-2 system. Such interpolation allows us to extrapolate data to a wide temperature range and obtain, firstly, information about the critical temperature from the intersection of the dependence χ − (1/T) and the calculated value χ cr (as shown schematically in Figure 1b). In accordance with this, the UCST for the PB-PS-1 system is 628 K, and for the PMMA-PS-2 system, it is 598 K. The calculated values for the concentrations at the critical point ϕ PS, cr by Equation (3) can be given as 0.64 and 0.91 for each of the systems, respectively. Second, extrapolation of the temperature dependence of the interaction parameter enables us to use data on the χ − (1/T) ratio to determine the corresponding compositions of coexisting phases with any temperature step. The phase equilibrium data for the PMMA-PS system and the temperature dependence of the interaction parameter are in good agreement with the data obtained by the turbidity point method for this system [24,25], taking into account the difference in the PMMA molecular weight.  Figure 6 shows the calculated data on the free energy of mixing and chemical potentials in the studied systems. It can be seen that the PB-PS-1 system is characterized by a more symmetrical concentration dependence for ∆ (Figure 6a). The extrema of the function is explicitly defined both for the case → 0 and for the case → 1. For the PMMA-PS-2 system (Figure 6c) for ~ 0.2, the first extremum of such a dependence is uniquely determined, and the localization of the second extremum at → 1 cannot be determined. When passing to the concentration dependences of chemical potentials (Figure 6b,d), according to the proposed method, in both cases, we can identify the selfintersection point and determine the corresponding compositions of the coexisting phases for the binodal curves.  Figure 6 shows the calculated data on the free energy of mixing and chemical potentials in the studied systems. It can be seen that the PB-PS-1 system is characterized by a more symmetrical concentration dependence for ∆G (Figure 6a). The extrema of the function is explicitly defined both for the case ϕ PS → 0 and for the case ϕ PS → 1 . For the PMMA-PS-2 system (Figure 6c) for ϕ PS ∼ 0.2, the first extremum of such a dependence is uniquely determined, and the localization of the second extremum at ϕ PS → 1 cannot be determined. When passing to the concentration dependences of chemical potentials (Figure 6b,d), according to the proposed method, in both cases, we can identify the selfintersection point and determine the corresponding compositions of the coexisting phases for the binodal curves.
For both systems, the concentration dependences of the iteration method were used to obtain data on the coexisting phase compositions for the binodal curves, and the coexisting phase compositions of the spinodal curves were calculated using Equation (10). This made it possible to construct generalized phase diagrams for both the PB-PS-1 and PMMA-PS-2 systems (Figures 7a and 7b, respectively). In both cases, the diagrams are characterized by UCST (gray dots), the position of which is localized in the region of high temperatures above the destruction temperatures of the components. However, compatibility information in this area, as well as in the area below the glass transition, is extremely important since blends in manufacturing processes can be obtained over a wide temperature range. In accordance with the difference in M w , the diagram for the PB-PS-1 system is more symmetrical than for the PMMA-PS-2 system. The calculated binodal (solid lines 1) and spinodal (dashed lines 2) curves correlate well with the initial experimental data on solubility (black dots) and accurately determine the position of homogeneous (I) and heterogeneous (II) regions, as well as the regions of metastable states (III) on the temperature-concentration field of the diagram. The dot-dashed line (3) shows the change in Alekseev's diameter (rectilinear diameter rule) when approaching the UCST. In both cases, this dependence is curvilinear, and for the PMMA-PS-1 system, this is more pronounced. At the same time, for the PMMA-PS-2 system, the right part of the diagram coincides to a lesser extent with the calculated values according to the binodal curve but coincides with the spinodal one. This can be explained by the error in the experimental determination of the solubility of the components in this region of the interferograms. For both systems, the concentration dependences of the iteration method were used to obtain data on the coexisting phase compositions for the binodal curves, and the coexisting phase compositions of the spinodal curves were calculated using Equation (10). This made it possible to construct generalized phase diagrams for both the PB-PS-1 and PMMA-PS-2 systems (Figure 7a and Figure 7b, respectively). In both cases, the diagrams the change in Alekseev's diameter (rectilinear diameter rule) when approaching the UCST. In both cases, this dependence is curvilinear, and for the PMMA-PS-1 system, this is more pronounced. At the same time, for the PMMA-PS-2 system, the right part of the diagram coincides to a lesser extent with the calculated values according to the binodal curve but coincides with the spinodal one. This can be explained by the error in the experimental determination of the solubility of the components in this region of the interferograms. It should also be noted that such an approach makes it possible to obtain generalized phase diagrams not only for systems with UCST but also for systems with a lower critical solution temperature (LCST). The technique was tested by us earlier on the classical system polystyrene-poly(vinyl methyl ether) (PVME) [22], which is characterized by the presence of LCST in the region of 370-470 K, depending on the PS molecular weight. In addition, the universality and simplicity of the mathematical apparatus of the Flory-Huggins theory and the presented approach make it possible to solve complex thermodynamic problems of simulation phase equilibrium diagrams in the presence of two critical solution temperatures. For example, for the PVME-water system [19], a generalized diagram was constructed with two LCST in the region of 290-310 K, and it was shown that part of this diagram represents the equilibrium between PVME and its complex with water, and the other part of the diagram is the equilibrium between PVME and water after the decomposition of the complex at an increase in temperature.

Conclusions
A method for calculating the coexisting phase compositions of binodal curves through the concentration dependences of chemical potentials in the framework of the Flory-Huggins theory was also proposed, using the example of the PS-PB and PS-PMMA systems. This made it possible to approach the construction of generalized phase-state diagrams with higher accuracy in cases where the critical point lies in the region of concentrated solutions due to a large difference in the molecular weights of the components. According to the proposed approach, there is no difficulty in determining the minima on the concentration dependence of the free energy of mixing since it is possible to fix the concentrations when the chemical potentials are equal.
The versatility and simplicity of the approach makes It possible to construct generalized phase diagrams containing binodal and spinodal curves for polymer-solvent, polymer-oligomer, polymer-polymer, and polymer-copolymer interactions characterized by both UCST and LCST.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.