Phase Behavior of Correlated Random Copolymers

In this work, we calculate Flory–Huggins phase diagrams for correlated random copolymers. We achieve it in two steps. At first, we derive a distribution function of two-letter A, B copolymer chains depending on the fraction of A-segments and AB-duplets. Then, we use the method of moments, which was developed by Sollich and Cates [Phys. Rev. Lett.80, 1998, 1365−1368] for polydisperse systems, to reduce the number of degrees of freedom of the computational problem and calculate phase diagrams. We explore how the location of transition points and composition of coexisting phases depends on the fractions of A-segments and AB-duplets in a sequence and the degree of polymerization. The proposed approach allows taking into account fractionation, which was shown to affect the appearance of the phase diagrams of statistical copolymers.


■ INTRODUCTION
Random copolymers, or statistical copolymers, are polymers composed of at least two monomer units connected in a more or less random manner. Random copolymers are ubiquitous in industrial applications. For example, poly(vinyl alcohol) (PVA) is used for the creation of water-soluble pouches and capsules. PVA is obtained by post polymerization modification from poly(vinyl acetate) (PVAc), in which the acetate groups can be either fully or partially transformed into alcohol groups. 1−4 Pure PVA films with a high degree of hydrolysis are brittle and difficult to dissolve because of a considerable degree of crystallinity. 5 Lowering the degree of hydrolysis or adding plasticizers leads to more flexible and soluble PVA films. To understand the conditions of stability of these films with respect to segregation and their physical properties, it is important to take into account a copolymer nature of PVA. 6,7 Other examples of statistical copolymers include carboxymethyl cellulose (CMC), 8 acrylonitrile butadiene styrene (ABS), styrene−butadiene rubber (SBR), membrane electolyte polymers like Nafion, 9 polyurethanes, 10 and many other commercial materials. 11,12 The ability to predict polymer material properties based on their chemistry, chain structure, and preparation method is crucial for product design and other industrial applications. The molecular structure of statistical AB copolymers is most commonly characterized by the composition (i.e., the fraction of A units). Composition totally describes an ensemble of sequences in the case when there are no correlations between the type of segments at different positions along the chain. In this case, copolymers are called random. If there is a correlation between the appearance of different types of segments at different positions along a chain, then statistical copolymers are called correlated. Correlated copolymers can be described macroscopically in terms of concentrations of duplets (AA, BB, AB, BA), triplets (AAA, AAB, etc.), etc. The more concentrations of n-tuplets are needed to fully describe the system, the less randomness there is in a copolymer sequence. Here, we consider the situation when only the fractions of Asegment and AB-duplets (which is equivalent to a fixed concentration of all duplets) are enough to fully characterize the copolymer. If the concentration of AB-duplets in a copolymer sequence is reduced, compared to a truly random copolymer, then such correlated copolymer is referred to as a blocky copolymer, meaning that the segments of one type tend to be arranged in blocks. In the opposite case when a copolymer sequence is enriched in AB-duplets, the copolymer is called alternating.
In the case of PVA, the degree of the blockiness of PVA depends on the synthetic route. 2,3 Saponification of PVA leads to blocky copolymers, while acetylation of previously hydrolyzed PVA produces more random sequences. The appearance of correlations during saponification is due to a larger reaction constant for hydrolyzing VAc segments whose neighbors have already been hydrolyzed. 3,13 Blocky or alternating copolymers are also common results of sequential polymerization. 14 In this case, blockiness arises when k AA ≫ k AB and k BB ≫ k BA , where k IJ is the reaction rate coefficient describing the addition of I-segment to the growing chain with segment J at the end.
Both models belong to the class of first-order Markov models and with the assumption that copolymerization is stationary (concentrations of monomers are kept fixed), both produce the same types of sequences.
The equilibrium phase behavior of statistical copolymers is rich and is not yet fully understood. Early works on the phase behavior of random copolymers 15−17 concentrated on considering Flory−Huggins mixtures or copolymers with different compositions (fractions of A-segments in AB-copolymer). As the Flory−Huggins parameter, χ, describing interactions between dissimilar segments, increased, an initially homogeneous mixture of copolymer chains separated into two phases with a different fraction of A-segments. The location of a spinodal point, with respect to separation into two phases, was predicted by Scott: 15 where ρ 2 and ρ 1 were correspondingly the second and the first moment of the distribution with respect to composition. Nesarikar et al. 17 went further and calculated phase diagrams with cloud points and fractions of A-segments in coexisting phases for short Bernoulli copolymers. They showed that as the Flory−Huggins parameter, χ, increased, separation into two, three, four, etc., phases occurred. Importantly, it was also shown that distributions with respect to the fraction of A-segments in coexisting phases, in general, had different shapes; in other words, fractionation took place. However, these predictions and calculations were made only for polymers with a small number of segments, N ≈ 10−30, and no correlations along the sequence.
Another approach to predict phase behavior of random copolymers was proposed by Shakhnovich and Gutin 18 and later developed by other authors. 19−25 It was based on Landau expansion of free energy of a copolymer melt in terms of an order parameter representing a deviation of a local composition from a global composition. Coefficients of expansion were calculated within a mean-field approximation and expressed through single-chain correlation functions averaged over all sequences. This approach was applied to correlated random copolymers by Fredrickson et al. 19 It was shown that initially homogeneous melt of blocky copolymers separated first into two macrophases upon an increase in the Flory−Huggins parameter. This transition was closely followed by remixing and forming one microphase separated phase with no long-range order. The period of the microphase had a strong dependence on the temperature and decreased as L ∼ (T s − T) −1/2 as temperature decreased (Flory−Huggins parameter increased). For alternating copolymers, a critical value of sequence correlation λ C was found, such that for λ < λ C a direct transition from the homogeneous state into microphase separated state was predicted, without macrophase separation. 19 Nesarikar et al. 17 pointed out that the one feature of the Shakhnovich−Fredrickson approach was that it did not take fractionation into account. It implicitly assumed that the distribution with respect to compositions (fraction of Asegments) had the same shape and only the mean value of composition had been changed. This assumption was too strong for nonsymmetric comopolymers. The further prediction was that for sufficiently short copolymers, N ≲ 60, the transition to microphase should be preceded by the coexistence of three macrophases.
Recently this discussion was continued by von der Heydt et al., 26 who considered a mixture of triblock copolymers, AAA, BBB, ABB, AAB, BAB, ABA, with overall composition f = 0.5 and varied volume fractions of different sequences to imitate Markovian sequence correlations. They showed that as the Flory−Huggins parameter, χ, increased coexistence of two Aand B-rich macrophases was followed by the coexistence of three phases one of which was the lamellar microphase. Microphase emerged as a shadow phase and was enriched in alternating sequences. So, both microphase separation and three-phase coexistence took place simultaneously.
The latest results show that it is important to take fractionation into account to make a correct prediction of phase behavior. 27 In this work, we aim at taking into account fractionation in the framework of the Flory−Huggins theory of blocky copolymers with realistic chain lengths. Therefore, we consider only the possibility of macrophase separation despite knowing about the existence and importance of microphase separation in statistical copolymers. This is done to get a solid reference point for a more refined picture including microphase separation, which may be developed in the future.
In this work, we use the method of moments proposed by Sollich et al. 28 for polydisperse systems. This method can effectively reduce the number of degrees of freedom of the polydisperse system; otherwise, the system consists of, an order of magnitude, 2 N different components and direct solution of the phase equilibrium equations is not possible. To use the method of moments, we derive the probability distribution function of copolymer chains with respect to the fractions of A-segments and AB-duplets. Then, we obtain Flory−Huggins phase diagrams of blocky copolymers and study the dependence of phase diagrams on the fraction of A-segments of the copolymer, chain length, and degree of correlations along the sequence (concentration of AB-duplets). At the end, we make a comparison of our results with the work of Nesarikar et al. 17 and Fredrickson et al. 19 The paper is organized as follows. First, we derive a distribution function for correlated copolymers. Then, we use this distribution to apply the moment method and obtain phase diagrams, volumes of coexisting phases, and their density distributions. We finish with a discussion of the results.

MARKOV COPOLYMERS
To derive a distribution function for Markov copolymers of the first order, let us first look at AB random binomial copolymers with chain length N, the number of A-monomers equal to N A , and their average fraction f = ⟨N A /N⟩ in the population of copolymer chains. The distribution function can be written as Let us consider now an infinite random AB sequence characterized by the fraction of A-segments, f. Then, we can write down an information rate that corresponds to entropy per monomer of such sequence 29,30 The information chemical potential of species A is then If there is a finite sequence of length N in equilibrium with this infinite system, then its grand potential depending on the number N A = σN of A-segments in this sequence is Also, the probability to have N A = σN segments is which is a binomial distribution. Now let us turn to the first-order Markov sequences. Information rate, in this case, is known to take the form 29 here, p A = f is the probability of a randomly chosen segment to be type A, p B = 1 − f is the probability of a randomly chosen segment to be type B, and n IJ is a concentration of IJ-duplets in the sequence. Concentrations of duplets satisfy conditions: n AA + n AB + n BA + n BB = 1, n AB = n BA = θ, n AA + n BA = f, and n BA + n BB = 1 − f, where f is the fraction of A units and θ is the concentration of AB-pairs. θ was referred previously as a block character by Moritani and Fujiwara 2 as it controls the average lengths of blocks when the composition of copolymer is fixed, ⟨l⟩ = (f)/ (θ). In completely random copolymer, θ = f(1 − f) and the average length of the A-block is then the sequence is depleted in AB and BA duplets and the copolymer is defined as blocky; the average lengths of A-blocks and B-blocks are larger than in random copolymer. Conversely, if θ > f(1 − f), the concentrations of AB and BA duplets are increased and the copolymer tends to be alternating, the average length of A and B-blocks is decreased. The maximum value the θ can take is min(f,1 − f).
An alternative parameter, which is often used to characterize the degree of correlations and the average block length in firstorder Markov copolymers, is a parameter λ, introduced by This expression depends on two parameters f and θ, which describe the concentration of A-segments and AB-pairs in the infinite sequence correspondingly. In analogy with the case of a random sequence, we calculate the chemical potentials of these pseudospecies from the information rate For a distribution function for a chain with the finite length N, the fraction of A-segments σ, and the fraction of AB-pairs t, which are in equilibrium with an infinite chain with the composition f and the fraction of AB-duplets θ, we get where Z is determined from the normalization condition If we reverse the Stirling approximation again in analogy with a binomial distribution, we can get (1 ) In Appendix B, we show another way to derive this distribution, which serves as additional support to the calculations presented above.

Macromolecules pubs.acs.org/Macromolecules Article
We would like to note that similar ideas were developed previously to characterize the local composition profile of compatible polymer blends in the course of macromolecular reactions and interdiffusion. 31,32 Example contour plots of the distributions are shown in Figure 1. We can see that at a fixed value of N and f, the width of the distribution projected on the σ axes decreases as θ increases. We expected to see this because the difference in compositions between different chains is expected to be larger when segments are arranged in blocks. If we compare two distributions with the same θ but different f, we can see that the distribution with f closer to 0.5 is broader, so variations in compositions are largest for the symmetric copolymer.

■ MOMENT FREE ENERGY
Let us now consider a melt of random copolymers with the distribution ρ(σ, t) and write down a moment free energy for it. 28 We start with Flory−Huggins free energy for the melt where V is the volume of the system, χ is a Flory−Huggins parameter describing interactions of segments of type A and B and Macromolecules pubs.acs.org/Macromolecules Article is the total volume fraction of A-monomers. R 0 (σ,t) is the distribution of chains with respect to the fractions of A-segments and AB-duplets in the parent phase. It is added here for convenience as far as it represents a term linear in density ρ(σ, t), which does not affect the phase behavior of a copolymer melt.
We also implicitly assume that volumes of A and B segments are equal to each other, both segments are flexible so that ρ(σ, t) is simply a volume fraction of the polymer with parameters σ, t. 33 Now we can apply the moment method 28 to reduce the number of degrees of freedom of the problem. We fix m first moments of the distribution ρ i = ∫ σ i ρ (σ, t)dσdt using Lagrange multipliers λ i along with the total volume fraction which is equal to 1 with Lagrange multiplier λ 0 and the constrained free energy is minimized with respect to the remaining degrees of freedom We note here that though we have two-dimensional (2D) distribution, depending on two parameters, ρ(σ, t), we fix only moments of the composition here. The reason for this is that the Flory−Huggins interaction term depends only on the first moment of the σ (fraction of A-segments), so fixing of other moments (for example, moments of t) does not affect the phase diagram. This was checked by direct calculations (fixing all moments up to the third order gives the same result as fixing only moments of σ). In case when the contribution of interactions into free energy depends additionally on θ 1 = ⟨t⟩, the average concentration of AB-duplets all moments of the distribution should be fixed (physically this case takes place, for example, when Flory−Huggins parameter describing interactions of two segments depends on the type of their neighbors, the model considered by Balazs et al. 34 ).
Minimizing F′ with respect to ρ(σ, t), we get In this approximation, all possible distribution functions in any coexisting phases belong to this family (eq 19). It is clear why we needed to include R 0 (σ, t) in the formula (eq 18), because the parent phase must be included in the family. The larger the number of fixed moments, m, the larger the set of functions for approximating distributions in daughter phases. Then, we substitute eq 19 into eq 18 and obtain an expression for the moment free energy (omitting terms which are linear in ρ i ) where χ′ = χN. This expression can be analyzed further as the free energy of a m-component system.

■ PHASE DIAGRAMS
In this section, we calculate phase diagrams using derived distribution (eq 15) and moment free energy (eq 20). We expect that upon an increase of the Flory−Huggins parameter, χ′, an initially uniform melt separates into two phases at a cloud point χ′ cloud , where a new phase with infinitesimal volume emerges. It is also expected that the spinodal point of the system is located at some value of χ′ s > χ′ cloud . Beyond the spinodal point, the a homogeneous state cannot exist as a metastable.
To construct a phase diagram and to determine the composition of coexisting phases from the moment free energy, we write down standard expressions for chemical potentials and osmotic pressures.
Conditions of phase equilibrium for coexisting phases α, β, γ, etc., are Additionally, there are conditions that amounts of ρ i are conserved where ρ i (0) is the ith moment of composition in the parent phase and the sum of all volume fractions of the coexisting phases equals to 1. These equations can be solved by the standard Newton method. For example, for three coexisting phases with m moments fixed, the set of unknowns is {λ 1 α ,λ 1 β ,λ 1 γ ,λ 2 ,..,λ m ,v α ,v β } and there are m+4 equations in total (see Appendix A for details).
To determine the first spinodal and cloud points, it is enough to fix only one moment. 28 For the spinodal point at which the initial homogeneous phase loses stability, analytic expression can be derived which is a well-known condition for the spinodal point in random copolymers. 15 To calculate the characteristics of coexisting phases above the cloud point, it is needed to fix more moments. However, it is not known in advance how many moments are needed to be fixed to produce a satisfactory approximation of the actual composition of the coexisting phases. We proceed by fixing an increasing number of moments at each step until the phase diagram does not change anymore for a given value of χN. Figure 2a shows the phase diagram of a copolymer with a fraction of A-units f = 0.3, fraction of AB-duplets θ = 0.09, and chain length N = 100. Figure 2b shows the dependence of the volume fractions of coexisting phases on χN. The binodal for separation of the initially homogeneous phase into two phases is located at χN = 56.2. One of these phases is the cloud phase with a composition equal to the initial composition of the system ρ 1 = 0.305 (it is denoted as α). It occupies nearly all of the system volume (see Figure 2b). Another phase is a shadow phase with a 2) corresponding to the instability point of the initial homogeneous phase. The decrease in the volume fraction of A-segments in the shadow phase is a consequence of the increase in its volume fraction. If the volume fraction of a phase is sufficiently high, its average fraction of Asegments is inevitably close to the volume fraction of Asegments in the initial phase because the initial distribution has one thin single peak. However, it is interesting to note that the compositional contrast between the phases above the spinodal, for example, at χN = 70 Δρ 1 ≈ 0.13 is larger than standard deviations of each phase (ρ 2 α − (ρ 1 α ) 2 ) 1/2 = 0.0829 and (ρ 2 β − (ρ 1 β ) 2 ) 1/2 = 0.0831 (noting that the standard deviation of the parent phase is a bit larger and is equal in this case to (ρ 2 0 − (ρ 1 0 ) 2 ) 1/2 = 0.0869) meaning that they can be distinguished. At χN ≈ 77, the second binodal with respect to the coexistence of three phases is located and a new shadow phase emerges (denoted as γ). The compositions of coexisting phases at this point are ρ 1 = 0.29, ρ 1 = 0.4, and ρ 1 = 0.64. At χN ≈ 86.4, the first shadow phase losses its stability, which corresponds to the second spinodal point (phase compositions are ρ 1 = 0.28, ρ 1 = 0.39, and ρ 1 = 0.54). Mathematically, this means that the determinant of the matrix composed of second derivatives of moment free energy with respect to fixed moments ρ i equals zero at this point in the shadow phase The next spinodal point is located at χN = 107.7, and here the shadow phase with the largest ρ 1 loses stability and separates into two phases. However, this fact is not reflected in Figure 2a. Also, the cloud point at which the fourth phase emerges is not shown. Below all diagrams are calculated with nine moments of compositions fixed. We discovered that to correctly predict the direction of change in the composition of a shadow phase above the first cloud point (for example, see Figure 2a) at least three moments should be fixed for copolymers with nonsymmetric composition f ≠ 0.5. To predict correctly, the location of the second spinodal point and volumes of coexisting phases, at least six moments need to be fixed. Also, to predict the location of the second cloud point, at least eight moments are required. Figure 3 shows the phase diagram of a copolymer with composition f = 0.3, blockiness θ = 0.15, and length N = 100, calculated with moment free energy depending on nine moments of the distribution. The increased value of θ, fraction of AB-duplets, compared to Figure 2a means that copolymer is less correlated. However, this value is still less than in a binomial copolymer with the same composition, i.e., θ = f(1 − f) = 0.21. One can see that the increase in concentration of AB-pairs leads both to an increase of the value of χN at which the first spinodal and the first cloud point are located. It is expected because an increase in θ leads to a decrease in the variance of the distribution, so increase of χN at spinodal point (eq 28). The same reason explains why the distance between consecutive spinodals on the diagram increases and the compositional  contrast decreases, though the composition of a shadow phase at the first cloud point does not change ( f sh = 1 − f). Figure 4 shows the effect of increasing N on the phase diagram, it is calculated for a copolymer characterized by parameters f = 0.3, θ = 0.09, and N = 300. Interestingly, we can see that as we increase N both the first spinodal point χ S and the first cloud point χ C slightly decrease and their difference Δχ = χ S − χ C essentially does not change.
Finally, Figure 5a shows a special case of a phase diagram for a system with a critical composition f = 0.5. It differs from the phase diagrams for asymmetric copolymers. The first spinodal point coincides with the cloud point, and the transition from one phase to two phases is continuous. The composition of the two coexisting phases is symmetric with respect to the line f = 0.5, and their volume fraction does not change until the third phase emerges (Figure 5b). The transition from two phases to three phases is discontinuous. As χN increases the volume fraction of a phase with the volume fraction of A-segments, f = 0.5 increases until the next spinodal line. The next transition is again continuous. These observations agree with the observations of Nesarikar et al. 17 for binomial copolymers. Figure 6 shows the dependence of cloud and spinodal points on the fraction of A-segments and chain length. As expected χ values at which spinodal and binodal points are located increase with the increasing asymmetry of the copolymer. A slight decrease of both χ S and χ C with N can also be observed. The spinodal converges to the limit derived by Fredrickson et al. 19 for N→∞ (see Figure 7) We also predict that binodal does not converge to spinodal in the limit N→∞ (Figure 8), and the difference between the Flory− Huggins parameter at the cloud point and the spinodal point, Δχ = χ S − χ C , is nearly independent on N. This is a new prediction. 19

■ DISCUSSION
In the previous section, Flory−Huggins phase diagrams for correlated copolymers are presented which were calculated using the method of moments. Here, we compare obtained results with existing works and propose future developments.
First, we want to note that our phase diagrams qualitatively look similar to those for short binomial copolymers considered by Nesarikar et al. 17 In that paper, phase diagrams were obtained by direct solution of phase equilibrium equations for all components in the system. Insofar as binomial copolymers are a special case of correlated copolymers, this demonstrates that the combination of distribution (eq 15) and the moment method produces consistent results.    Our prediction for a spinodal converges to the expression derived by Fredrickson et al. 19 for infinitely long chains. We also showed that the distance between binodal and spinodal is nearly independent of N in contrast with previous predictions. 19 We predict that the compositional contrast between coexisting phases is larger than the width of the distributions of these coexisting phases, so in the framework of Flory− Huggins theory, phases can always be distinguished. It would be interesting to verify our predictions with respect to contrast in the composition of macrophases in simulations (at least for short sequences), for example, using the method proposed by Houdayer and Muller 35 extended to simulate asymmetric copolymers.
In the present work, we consider the simplest possible case of the copolymer in a melt with a fixed chain length. However, the method of moments allows taking polydispersity in chain length into account as well. The only thing which is needed for it is a distribution function depending additionally on the chain length, ρ(σ, t, N). In the case of the PVA−PVAc system synthesized by the postmodification of PVAc, this distribution function is obtained by the product of distribution function describing polydispersity of PVAc and the distribution function with respect to the fractions A-segments and AB-duplets (eq 15). Additionally, for the PVA−PVAc system, it should be taken into account that the volume of the VAc monomer unit is 2 times larger than the volume of the VA monomer unit, which leads to polydispersity in length even in the case when the initial PVAc is perfectly monodisperse. This effect is also easy to account for using the presented approach. Both polydispersity in length and segment asymmetry increase incompatibility. Interestingly, our preliminary calculations show that the segment asymmetry has a larger effect on the phase diagram than the chain length polydispersity (for Poisson distribution of initial PVAc).
The proposed approach can also be generalized to predict the phase behavior of mixtures of Markov copolymers with a plasticizer (solvent), which is an important industrial problem. It is well known that the phase behavior of these mixtures is strongly affected by the polydisperse nature of these materials, especially for blocky copolymers. 28,36−38 To derive a distribution function for the first-order Markov copolymers, we used an approach related to the large deviation theory, 30 which allows us to do it easily and additionally go beyond a Gaussian approximation. 39 To understand the effect of large deviations in the initial distribution on phase behavior in the framework of our model, we compare (see Figure 9) phase diagrams obtained with the method of moments for distribution (eq 15) (black curves) and the Gaussian approximation of this distribution (red curves) in a case f = 0.3, θ = 0.09, and N = 100. As far as the variance for both distributions is the same, the first spinodal point is the same in both cases, it is shown by a dotted line in Figure 9. The composition of the shadow phase at the cloud point is also the same in both cases and equals 1 − f. The location of the first cloud point is different (shown by dashed black and red lines). In the case of Gaussian distribution, it is located just below the spinodal in contrast to the case of parent distribution (eq 15) considered in this paper. The contrast in composition between cloud (α) and shadow (β) phases is small in the case of Gaussian distribution. It is smaller than the width of the distribution within each coexisting phase, so such coexisting phases cannot be distinguished. In contrast, for the case of the distribution (eq 15), phases can be distinguished as we discussed above. This means that difference here is not just quantitative but also qualitative. The last difference is in the location of the second spinodal point. In the case of Gaussian parent distribution, the first shadow phase loses stability at χN ≈ 73.3 (second spinodal point shown with a red dotted line in Figure 9), which is a much smaller value of χN than in the case of parent distribution considered in this paper. We did not go beyond the second spinodal for the parent Gaussian distribution because we were unable to find a converging solution for three-   phase coexistence in this case. Concluding, we can suggest that Gaussian approximation in this particular case gives qualitatively different predictions for phase behavior.
In the present work, we use the Flory−Huggins model to describe a copolymer melt. There are several limitations to this model. One of them is that it does not allow to consider the possibility of microphase separation. It is expected that the picture including microphase separation is the following. Upon increase in Flory−Huggins parameter initially, homogeneous melt separates into two coexisting macrophases. 19 At further increase in Flory−Huggins parameter, first, microphase separated phase coexisting with two macrophases is expected to emerge, 26 then it is expected that macrophases will be absorbed by a microphase. The type of microphase was predicted to depend on the composition of the copolymer 21,25,40 and the degree of incompatibility. The possibility of coexistence between different ordered structures was also predicted. 25,41 All of these results were obtained with the assumption that distribution in all coexisting phases had the same shape and only the average composition was different. It would be interesting to see in the future theory developments accounting both for microphase separation and fractionation for the case of asymmetric random copolymers.

■ CONCLUSIONS
In this paper, we derive a probability distribution function for a blocky copolymer, with first-order Markov correlations along the sequence. Then, we used this distribution with the method of moments by Sollich et al. 28 and Flory−Huggins theory to obtain the phase diagrams of blocky copolymers.
This combination of methods allows one to calculate Flory− Huggins phase diagrams for copolymers characterized by arbitrary length and blockiness. Qualitatively obtained phase diagrams look similar to those for short binomial copolymers. 17 As the final message, we would like to emphasize the importance of characterizing Markov copolymers (those obtained in stationary conditions) by both fractions of Asegments and AB-duplets. This work shows that the phase behavior of two copolymers with the same composition but a different fraction of AB-duplets can be very different. So, characterization of the structure of copolymers, such as PVA, in practical settings should include the determination of this parameter.

■ APPENDIX A
In this section, we include computation details. Consider the case of moment free energy with m moments fixed and two phases in equilibrium. Then, to determine the compositions of these phases, one needs to solve a system of equations , ..., It is solved numerically with the Newton method i k j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z All calculations were made using scripts written in Python with the mpmath package allowing arbitrary precision mathematics. The need in arbitrary precision mathematics arises because the matrix (eq 32) in general is poorly defined. All integrals were calculated with the trapezoidal rule and uniform discretization along the composition axis. For factorials, Stirling approximation including a square root was used π ! ≈ · n n n e 2 ( / ) n . The step along χN axes was varied to ensure convergence of the Newton solver and maximization of the computation speed.
The transition from discrete distribution to continuous, as well as following discretization, introduces errors upon calculation of integrals. These errors lead to relatively large dependencies of λ i parameters on the degree of discretization (of order 1). These differences increase with χN and decrease with an increasing degree of discretization. Differences in λ i , however, produce much smaller differences in free energies calculated with different discretizations, which are of order 10 −5 for discretizations up to 500. However, the difference between the energies of two and three phases is of order 10 −7 when all calculations are done with the same discretization. As a result, the free energy of two-phase coexistence calculated with larger discretization may be smaller than the free energy of three-phase coexistence calculated with smaller discretization, even at the second spinodal point when the first shadow phase losses stability. So, it is recommended to keep discretization fixed during calculations. Values of all transition points depend to some extent on discretization. However, the difference between discretization 100 and 500 is of order δχN ≈ 0.1 and differences tend to decrease with an increase in discretization. As far as computational cost is concerned, this increases strongly with an increase of discretization (∼N discr 2 ). We used discretization 100 for calculations of all phase diagrams.

■ APPENDIX B
Consider again the binomial distribution In this expression, f N A (1 − f) N−N A gives the probability of a specific sequence ABABB... with total length N and the number of A-segments equals to N A . Also, is the total number of sequences which has the same description, i.e., have a total length N and the number of A-segments equals to N A . So, if we look at the distribution for a blocky copolymer we can see that it has the same structure as the binomial distribution and the probability of a specific sequence, ABABB..., in which we calculated both the number of A-segments, N A , and the number of AB-duplets, N AB equals to (1 ) So, if we calculate the probability of ABBB sequence (it is short here for simplicity), we have A AB BB BB 2 2 (38) where P IJ are the corresponding elements of the transfer matrix (eq 32). From the formula above (eq 37), we obtain So we have a discrepancy between p(ABBB) and p* (ABBB). The contradiction is resolved if we take into account that in eq 37 no start of the sequence was ever specified and it must be "cyclic" without cyclic symmetry. So to obtain the probability for a linear sequence, we need to specify a starting segment and remove one duplet adjacent to it (in this particular case the BAduplet) = * · p p p P (ABBB) (ABBB) A BA (40) However, these end-chain effects will be unimportant in the limit of large N, N A , and N AB . The second part of the distribution, the number of sequences with the same description N, N A , and N AB , can be calculated in the following way. Assuming that we have a string of A-s with length N A and a string of B-s with length N − N A . Then, if we want to mix them in a way such that there are precisely N BA + N AB = 2N AB boundaries between A and B-blocks, we need to find a number of ways in which one can split a continuous A-string into N AB blocks and for each of these splits find a number of ways the B-string can be split into N AB blocks. This means that the total number of combinations is This expression coincides with the number of sequences with the same description from the distribution in eq 36. In this case, end-chain effects are also neglected as the number of boundaries 2N AB is even, implying that all beads are located on a "circle".
It is interesting to note that distribution in eq 36 converges exactly to binomial distribution in the case when θ = f(1 − f) and the summation over N AB is taken. This happens because for binomial distribution there is no difference between a linear chain model with ends and the circle model considered above without ends, as there is no correlation between the probability of appearance of different segments.