Theoretically predicting the solubility of polydisperse polymers using Flory–Huggins theory

Polydispersity affects physical properties of polymeric materials, such as solubility in solvents. Most biobased, synthetic, recycled, mixed, copolymerized, and self-assembled polymers vary in size and chemical structure. Using solvent fractionation, this variety in molecular features can be reduced and a selection of the sizes and molecular features of the polymers can be made. The significant chemical and physical dispersity of these polymers, however, complicates theoretical solubility predictions. A theoretical description of the fractionation process can guide experiments and material design. During solvent fractioning of polymers, a part of the polydisperse distribution of the polymers dissolves. To describe this process, this paper presents a theoretical tool using Flory–Huggins theory combined with molecular mass distributions and distributions in the number of functional groups. This paper quantifies how chemical and physical polydispersity of polymers affects their solubility. Comparison of theoretical predictions with experimental measurements of lignin in a mixture of solvents shows that multiple molecular features can be described well using a single set of parameters, giving a tool to theoretically predict the selective solubility of polymers.


Introduction
The worldwide production and use of polymeric materials increases, and without intervention, is predicted to double by 2030 [1] and almost quadruple by 2050 [2].Until the end of 2050, the total global cumulative production since 1950 is expected to increase to, for example, 26 Gt of resin, 6 Gt of polyphtalamide fibers, and 2 Gt of additives [3].In that year, 20% of the total oil consumption is expected to be used by the plastics sector, [2,4] and 15% of the global annual carbon budget [2].The use of fossil fuels to produce these materials is not sustainable and has an undesirable contribution to global CO 2 emissions.
At the same time, awareness of the environmental impact of plastics is growing, as reflected by an increasing number of scientific publications [4] and media coverage [5].Part of the solution to the problems with fossil-fuel-based polymeric materials is the development of new bio-based polymeric materials.Carbon dioxide is captured in bio-based feedstock, [6] and using bio-based polymers instead of fossil feedstock-based polymers reduces greenhouse gas emissions [7].However, despite the attention for the environmental impact of plastics [4,5] and the high substitution potential for bio-based alternatives, [8,9] polymeric materials are still largely based on fossil-feedstock [3].
Improving the understanding and tunability of bio-based materials and the processing of bio-based materials can help implement these materials by providing options offering more desirable properties while maintaining the environmental advantages.A challenge for many of these polymers is that they vary in molecular features [10].One example is the bio-based material lignin.Lignin is a structural material consisting of aromatic subunits present in plants, that vary in microscopic composition depending on the source and used processing methods [11].Reviews show the advantages of potential applications of lignin, [12][13][14] including fibers, [15][16][17][18] composites, [19][20][21] and nanoparticles [22][23][24].Despite its potential, lignin is still hardly used to manufacture materials [12,14,25].Reasons for this are lignin's varying composition, structure, and polydispersity, and as a result variable and mostly unpredictable material properties [26].
Like lignin-based materials, other applied polymeric materials are typically polydisperse.For example, the degree of polydispersity and the shape of the molecular mass distributions of polystyrene affect its melt rheology.[27][28][29] Also, widening the polydispersity improves the processability of polypropylene due to its effect on shear thinning [30,31].Another effect of polydispersity is its influence on the glass transition temperature of poly(methyl acrylate) [32] and polystyrene [29].Moreover, dispersity also affects the properties of copolymers, such as the order-disorder transitions of block copolymers [33,34] and self-assembly in bulk or solution [35].These effects of polydispersity, together with the effects on many other properties, indirectly affect the many applications of polymers in, for example, materials, coatings, food, and health care.Each of these applications could benefit from more well-defined molecular features.
There are multiple ways to control the polydispersity of polymers.Copolymer polydispersity can be controlled using catalysts, using reagents, or by regulation during the polymerization [35].Polydispersity can be increased by combining polymers with different distributions [33,35].Polydispersity can be decreased by selectively obtaining a subset of polymers using the variation in their properties using for instance membrane technology or chromatography.One can also take advantage of the differences in solubility, mobility, melting point, or weight to fractionate polymers [36,37].Solvent fractionation is a technique used to optimize the production and characterization of materials where a subset of polymers is selectively dissolved.The wide variety of possible solvent combinations makes fractionation using solvents tunable and specific.Solvent fractionation is beneficial for both analytical and preparative work, allows for large-scale fractionation, and has a low cost compared to many other separation methods [36,37].
Various publications [36,[54][55][56][57][58][59][60] discuss possible methods to account for variation in polymer length in other polymeric systems.Accounting for polydispersity directly in numerical computations, using, for example, successive substitution is sometimes possible, but often computations fail to converge [61].Alternatively, a continuous mass distribution can be simplified to a finite number of lengths.However, computation becomes unmanageable when a large number of lengths is used and introduces systematic errors for a lower number of lengths [56].A third option is to use perturbation methods, but these can only be used if the width of the chain-length distribution is small [62].Although each of these methods is useful for some instances, they do not enable the calculations for the known high physical and chemical polydispersities of solubilized lignin.
Here, we study solvent fractionation using Flory-Huggins theory [63,64].Huggins already compared the solubility of polymers with different chain lengths in 1967 [65].In this comparison, he, however, did not account for the interplay between the different-sized polymers in polydisperse samples.Cantow [36] and Krigbaum and Geymer [54] use extra correction terms in their expressions for the chemical potentials to account for polydispersity.Although these chemical potentials with corrections can fit experimental data, they cannot be used to predict solvency or explain its theoretical background.Chemical variation in nonuniform and chemically heterogeneous copolymers have been considered by averaging, [55,58,59] or assuming no interactions between polymers [60].Furthermore, each of these methods do not account for the fractionation effect in which smaller polymers dissolving at higher concentrations and larger polymers dissolving at lower concentrations affect the solubility of polymers of a certain chain length.
Recently, we demonstrated theoretical predictions agree with measurements of the yield of solubilized lignin as a function of the composition of a solvent mixture of methanol and ethyl acetate [66].The second application in the appendix of this publication provides extra support for the inclusion of physical and chemical polydispersity to predict the solubility of a different lignin type in a different solvent mixture.This paper presents the theoretical method used for these theoretical predictions.Here, we show how physical and chemical polydispersity affects solubility and how this solubility affects the molecular features of the dissolved polymer.Also, we provide detailed theoretical derivations and discuss their accuracy.

Theory
The Hansen solubility parameters [67] and Flory-Huggins interaction parameters [63,64] are measures for the interaction between polymer segments and solvent molecules.In this paper, we relate these parameters to compare theoretical and experimental yield of solubilized polymers.We define size dispersity as the variation in the chain length of the polymers and chemical dispersity as the variation in the amount and types of functional groups.The yield is here the percentage of polymers that dissolve relative to the total amount of polymer.
The Hansen solubility parameters compare the interaction energies of polymers and solvents to estimate miscibility.According to Hansen solubility theory, a polymer and solvent mix when the Hansen solubility parameters of these components are similar enough.Hansen solubility parameters can be combined to calculate a distance (R a ) in Hansen space [67].For a mixture of polymer and solvent: In this equation, δ dp , δ pp and δ hp and δ ds , δ ps and δ hs are the dispersion, polarity and hydrogen-bonding energy of the polymer and the solvent respectively.The relative energy difference (RED) is defined as: When the Hansen solubility parameters are similar, the location of the mixed components in Hansen space is similar and the components mix.For infinitely long, monodisperse polymers in contact with a solvent, the interaction radius (R 0 ) denotes the maximal distance in Hansen space where the components are still miscible.Using that these polymers dissolve if the RED ≲ 1 and can phase separate if the RED ≳ 1; R 0 can be determined using the solubility in different solvents.The considered experimental system is lignin dissolved in a combination of methanol and ethyl acetate.The used Hansen solubility parameters of methanol and ethyl acetate are respectively δ dm = 15.1, δ pm = 12.3, and δ hm = 22.3 and δ de = 15.8, δ pe = 5.3, and δ he = 7.2 [67].The Hansen solubility parameters of the mixture δ is can be calculated from the volume average of the individual Hansen solubility parameters δ im and δ ie [67].For lignin, we used δ dl = 21.42,δ pl = 8.57 and δ hl = 21.80 and R 0 = 13.56 in accordance with Novo et al [49].Other reported values for the Hansen solubility parameters of lignin are close to these values [50,51].
Within Flory-Huggins theory, [63,64] the parameter χ is a measure of the interaction energy between molecular entities.There is a strong parallel between the Hansen solubility parameters and the Flory-Huggins theory parameters for infinitely long monodisperse polymers.To calculate the RED from χ, we use: with χ c the χ at the critical point.For infinitely long, monodisperse polymers in a solvent, a similar logic as for the Hansen solubility parameters can be used [68].These components mix at all concentrations if χ < χ c and phase-separate at all concentrations for χ > χ c .Equations (3) follows directly from the quadratic dependency χ ∝ R 2 (and χ c ∝ R 2 0 ).[69,70] Variations in the dependency of χ on δ slightly modify equation (3), for example when χ is split in an enthalpic and entropic part [71,72].Using such variations of equation (3) in the theory presented below can for example be useful to account for temperature changes, but we do not consider this here for simplicity.
Polymers that are not infinitely long and monodisperse, can have a RED-range with varying solubility.Flory-Huggins theory has the advantage that this dependency can be predicted based on thermodynamic properties.The maximal miscibility or solubility of polymers with a finite size can be described by the binodal concentrations predicted using Flory-Huggins theory.Using the relation between the RED and χ in equation ( 3), full miscibility is achieved when the RED< 1, whereas the solubility is limited for RED⩾ 1.Also polydispersity affects the miscibility in multiple ways.The amount of variation in the polymer chain lengths and the number of functional groups affect the mixing entropy and enthalpy, which affect the miscibility.Furthermore, R a and χ c depend on the size and number of functional groups in the polymers.For calculations including polydispersity, it can be asserted that components are fully miscible for RED < 1, by using the least soluble polymers to determine R 0 and χ c .Since the least soluble polymer dissolves, it is certain that all other polymers also dissolve.The least soluble polymer can be the polymer with the highest length or the highest interaction energy of the present chemical entities.In this paper, we use two ways to determine R 0 and χ c .In sections 3.1-3.3,we calculate R 0 and χ c in all computations using the limit of infinitely long monodisperse polymers, so polydisperse and monodisperse results can be compared.In section 3.4, we use as a reference that the yield is 70% at RED = 1, such that the theoretical and experimental results can be compared.
The equilibrium concentrations in coexisting phases α and β, which are in thermal and mechanical equilibrium, can be calculated for each component k using that: In each of the phases the chemical potentials following from Flory-Huggins theory are given by: [63, 64] Here, M i is the chain length of component i, χ ij is the interaction parameter and is the volume fraction of components of type i in phase α.
In Flory-Huggins theory, the chains with different chemical compositions and lengths can be considered to be different polymers.With many molecular variations, in principle, many coexisting phases can coexist.Due to the chemical similarity of the polymers, we assume here that there are only two coexisting phases, a solvent-rich phase S, with a low volume fraction of polymer (ϕ S p ) and a polymer-rich phase P, with a high concentration of polymer (ϕ P p ).In the case of a (polydisperse) polymer and a solvent in two phases P and S, thermodynamic equilibrium therefore requires: as follows from equations ( 4) and ( 5).In these expressions, i, j and k are the indices of the different molecular variations.Because the number of phases is lower than the number of components, the concentrations are not only determined by the chemical potentials but also by the amount of each component present.These relations are used in the next section to derive expressions for the molecular features of dissolved, polydisperse polymers.In mean-field theories, the effective interactions of polymers with homogeneously, mean-field mixed polymers or solvents with different interaction parameters can be determined by taking an effective average of surrounding polymeric groups [59].In reality, this is, however, often not straightforward.
The number and location of functional groups vary for each polymer, and the interaction energies also depend on the interaction energies of neighboring groups.For example in systems with random copolymers, these parameters will never be exact.For the theoretical description of the solubility of lignin, we chose to limit the number of parameters by assuming that there are pairs of different monomers with a set interaction energy χ f and indistinguishable monomers for which χ = 0.The total interaction parameter χ of a polymer is the result of the mean-field averaging of the interactions of the different groups.For example, between i and j: with θ i the relative amount of monomers with functional groups in polymer i and θ j the relative amount of monomers with functional groups in polymer j.Interactions between differing polymers increase χ ij , while the interactions between identical polymers decrease χ ij .The same logic is used for the interactions with the solvent.Although neglecting the differences for the different χ's of different functional groups neglects some of the system's complexity, it allows for the inclusion of the other often neglected system characteristics of chain length and chemical polydispersity, and reduces the number of difficult-to-determine parameters.
By considering variations in the chain size and number of functional groups in polymers as different polymer types i, polymer chain length and distributions in functional groups can be considered in the calculation of the equilibrium concentrations in the polymer-and solvent-rich phases.The resulting binodal concentrations can be transformed to the yield by dividing the total amount of dissolved polymer with the total amount of polymer.In the next sections, we use this concept to derive equations for the yield as a function of the RED and apply those equations to quantify lignin dissolution in solvent-mixtures.

Results and discussion
Many polymers are disperse in chain length and chemical composition [10].This polydispersity can be included directly in numerical computations.However, solving these equations becomes extremely challenging for wide chain-length and chemical distributions.Sections 3.1-3.3discuss our general approach of including the size and chemical polydispersity in the calculations of the equilibrium yield of a complex mixture of polymers in a solvent.The derivation of these yield equations consists of three steps.First, we derive approximate analytical expressions for the equilibrium concentrations of the different species of polymers in the coexisting phases, giving upper-and lower-bounds for these concentrations.Two limiting results are derived.In section 3.2, these individual concentrations are combined and used to calculate the total dissolution yield.Finally, examples are given of distributions in polymer chain length and functional groups and their effect on the dissolved fractions in section 3.3.The theory derived in these sections is applied to lignin in section 3.4.

Solubility limits of polydisperse polymers
The first step in finding the yield as a function of the RED based on Flory-Huggins theory is finding the equilibrium concentrations of the coexisting phases.When a polydisperse polymer partially dissolves in a solvent, two phases will form: one phase with mostly polymer compounds and some solvent and another phase with mostly solvent and a bit of polymeric material.To find the relative amount of polymer that dissolves, we need the polymer concentrations in both phases.When a mixture consists of a limited number of components, binodal concentrations can be found directly numerically.However, when there are more components in a system, numerical determination of the equilibrium concentrations becomes increasingly challenging due to the high number of variables and equations that must be satisfied.Polydisperse samples have the same problem because polymers with different chain lengths and functional groups have to be treated as distinct polymers.To solve this problem, we determine the equilibrium concentrations by combining assumptions in two limits in which the equations can be solved, giving a range of possible values for the concentrations.
In the considered system, a small amount of polymer is added to a large amount of solvent (ϕ s ≫ ϕ p ).Even if all polymer dissolves, the concentration ϕ S p of polymer in the solvent is small.The range of possible concentrations in equilibrium of polymer in the solvent is 0 < ϕ S p < ϕ p .This can be used to determine the range of possible binodal concentrations of each polymer variation as a function of the RED.First, the maximal concentration as dictated by the considered system is combined with the chemical potentials from Flory-Huggins theory to estimate ϕ P s .In this limit, the polymers that are on the verge of being or not being dissolved can be described with M max and χ max ps ; the chain length and interaction energy with the solvent of the least soluble polymers.So the polymer-rich phase in this limit consists of only two components.As long as the total amount of polymer is small compared to the amount of solvent, ϕ P s can be approximated with a expansion of the logarithmic-term in equation (6b) around the maximal concentration.Our earlier publication [73] gives additional information on the validity and derivations of this approximation.Using the upper limit of ϕ S p , M max , χ max ps , ϕ p and a single, numerical computation of equation ( 4) in equation ( 6) gives: Figure 1.The maximal expected solvent concentration ϕ P s max obtained from equation ( 9) and computed numerically for a two component system as a function of χps (a) and Mp (b).The plot of equation ( 9) uses M max p = 10 4 and ϕp = 0.03.The chain length of polymers doubles with every next numerical curve in (a), while χps is varied in (b).On the right-side of figure (b), the red lines indicate the limit used in the equations, which are correct for higher chain lengths.For smaller Mp, concentrations approach the critical concentration (purple triple-dotted-dashed curve), which is the highest possible concentration for small Mp.
Next, eliminating ϕ S p from the equations and solving for ϕ P s provides the expression for ϕ P s : Figure 1 shows how the dependency for the maximal solvent concentration in the polymer-rich phase (equation ( 9)) compares to numerical computation for multiple values of χ ps and M for a two-component system.As expected, as long as polymers are not too small, ϕ P s is close to the numerical values (right-side of figure 1(b)).Closer to the critical point (lower χ max ps or M p ), the numerical values can be higher, but this is anyway outside of the considered range of ϕ S p values.As long as M max is chosen high enough, the effect of only using one value for M max p is small, even when the resulting concentrations are used for smaller M p .The maximal value of ϕ P s determined above can be used to determine ϕ S p .Again, equation ( 6) can be simplified to an expression directly giving the ratio of the concentrations of each of the polymer variants k.For systems with large polydispersity, each individual polymer concentration ϕ P k is small in comparison to the total amount of polymer.For small ϕ P k , the summations ∑ m̸ =k can be replaced with ∑ m without a large effect on the calculations, due to the relatively small effect of the small concentration of polymer type k.This makes it possible to approximate equation (6b) as: The summations are now independent of k, making it possible to simplify equation ( 10) further.For the effective masses, we used the approximation of Enders, [55] Vp the relative volume of polymer i in proportion to the total volume of all polymer.The effective interaction energies between polymer k and the other polymers is assumed to be ⟨χ kp ⟩ = ∑ i χ ki Θ i and the effective interaction energies between all polymers is assumed ⟨χ Using these approximations for the summations in equation (10) gives: ] .
(11) Equation ( 7) defines χ ki and χ ij .The advantage of this approximation is that the right-hand side of the equation is now independent of the concentrations ϕ S k and ϕ P k .As a result, the ratio of the concentrations ϕ S k /ϕ P k can be determined directly from the equations, for a given ϕ P s and ϕ S p .The problem here is that these concentrations are yet unknown.We do, however, have a range of possible concentrations.Above, we showed that the maximal concentration of polymer in the solvent-rich phase is ϕ and the maximal concentration of solvent in the polymer-rich phase ϕ P s is given by equation ( 9).Concentrations cannot be negative, so the minimal concentrations are ϕ S p min ≡ 0 and ϕ P s min ≡ 0. Using these limits in the calculations, gives two limiting dependencies for the concentrations as a function of the RED.
Combining equations ( 9) and ( 11) and the effective parameters results in an analytically solvable range of possible ratios ϕ S k /ϕ P k .These concentrations are used as input in the next section to derive the equilibrium polymer solubility.

Calculating the yield of solubilized polymer
When polymer is added to a solvent, they mix completely or form two or more coexisting phases.In section 3.1, we derived expressions for the equilibrium volume fractions of the polymer in the polymer-rich phase and in the solvent-rich phase.Here we convert these concentrations to the yield; the fraction of polymer that dissolves relatively to the total amount of polymer.This yield can be compared directly to the experimentally measured yield.First we derive the dissolving fraction η k of a single polymer variation k.Using equation (11) to eliminate ϕ P k in V k gives: with V P and V S the volumes of the polymer-rich and solvent-rich phase, which depend on the yield η, which is yet to be determined.With the volume fraction ϕ S k now known, the fraction dissolved of a specific species η k can be determined: When the volume of polymer added to a solvent is low compared to the solvent volume, the volume fraction in the solvent is low even when all the polymer dissolves.As a result, the range of yields where the approximation is accurate can be larger than the range of concentrations where the approximation is accurate.Figure 2 gives some examples of this effect.This figure shows that the smaller polymers dissolve completely before reaching a concentration close to the critical concentration.For longer polymers, the critical concentration is lower than the assumed maximal concentration, so this maximal concentration is never reached and the approximation is correct for all concentrations and all yields.following from the approximations in the maximal and minimal limit as discussed in section 3.1.The curve for Mp = ∞ is the dependency obtained directly from the Hansen solubility parameters (η = 1 for RED < 1 and η = 0 for RED > 1), which is valid in the limit of infinitely long polymers.The gray curves represent the numerical results of the Flory-Huggins theory equations ( 4) and ( 5).
The fraction of polymer dissolving of all types (yield) can be found by summation over all the different species: with Vp .The amount of polymer relative to the amount of solvent in the system is small, so V P and ϕ P s are small.Thus we can assume V P V S is not affected by solvent entering the polymer-rich phase, or V P ≈ (1 − η) V p and V S ≈ V s + ηV p , giving: If dk is small enough, this can be converted to an integral: The analytical limits for the concentrations can be used in equation (15) or equation ( 16) to obtain a solvable relation for η.Whether the solution can be solved analytically or only numerically depends on how Θ k , M k , χ sk and χ km depend on k.To find the yield of solubilized polymer, we solved equation ( 15) numerically.

Distributions in polymer chain lengths and functional groups
The final step in the calculation of the polymer yield is the inclusion of the variation in chemical groups and chain lengths of the polymers.We demonstrate this procedure with a couple of examples.As a first example, for a monodisperse system with polymers with an interaction energy corresponding to a functional group concentration of θ 1 (resulting in the same χ ps for all polymers) and chain length M 1 : In this case, the molecular mass M k and number of functional groups per segment θ k are independent of the index of the polymer k.Insertion of Θ k = 1, χ ps and M 1 in to equation (11) gives the expected dependency for monodisperse polymers in a solvent in the limit of small concentrations.Using equation (17) in equation ( 16) gives the yield for monodisperse polymers dissolved in a solvent.Secondly, we apply a more realistic log-normal distribution in chain length: Constant a converts the measured molecular mass to the chain length of theoretical, dimensionless beads.The quantities σ and µ are constants directly related to the molecular mass distribution of the polymers, which can, for example, be obtained using a fit through such a distribution.Figure 3(a) gives an example of the effect of chain length polydispersity.Here, we used an average polymer size ⟨M⟩ of 100 monomers and varied σ.Smaller polymers have a higher entropy per segment than longer polymers, resulting in a higher yield at higher RED for more polydisperse solutions.At the same time, longer polymers have a lower entropic contribution per segment on the other side of the distribution.They will only dissolve for lower RED than polymers with the average chain length.Similarly, variations in the number of functional groups per segment can be included.As an example, we will consider here dissolving of polymers for which there is a normal distribution in the number of functional groups per segment: Each functional group increases χ f with ∆χ = 2.The average number of functional groups per segment is in these calculations θ = 1 3 and the width of the distribution around this value is varied.The variables σ and µ can be obtained from a fit through distribution in functional groups from experimental measurements.Figure 3(b) gives examples of the effect of a distribution in the number of functional groups per segment.It should be noted that for each of the distributions, the RED is normalized using the monodisperse distribution (and M = ∞) to be able to compare these graphs.If the least soluble polymer had been used for the normalization in each calculation, η would reach 100% at RED = 1 for all σ, but the RED would correspond to different solvents at this value.
The results show that the interactions of some of the polymers result in a total fraction of polymer that dissolves that is higher at higher RED, while lower at lower RED.The polymers with more attractive interaction energies than the average interaction energy will dissolve already at a higher RED, while the ones with more repulsive interactions with the average interaction energy dissolve at lower RED than their more monodisperse equivalents.Comparing the effect of physical and chemical polydispersity shows that both increase η for high RED, due to the smaller and less repulsive polymers.A difference between the two effects on the η/RED-dependency is that chemical group polydispersity reduces η more for lower RED, due to the larger effect of higher χ k on the enthalpy than the effect of higher M k on entropy for polymers that have a high average chain length.The fractionation of physical polydispersity in figure 4(a) could for example be applied to reduce the width of the distributions in molecular mass of the earlier mentioned polystyrene, polypropylene or poly(methyl acrylate), [27][28][29][30][31][32] or other physical polydisperse polymers.The fractionation of chemical polydispersity in figure 4(b) could be applied in the preparation of the earlier mentioned copolymers [33,34] or polymers for self-assembly, [35] or other chemically polydisperse polymers.
The different predictions of how η dependents on M k and θ k can be used to learn about the properties of the dissolved and undissolved polymer fractions.Using η in equation (15a) for η k , makes it possible to directly calculate the fractions of dissolved polymer for each polymer variation.Figure 4 shows the resulting distributions in each of the phases of the log-normal distribution in chain lengths (a) and normal distribution in the number of functional groups per segment (b).For example the steepness of the boundary between dissolved and undissolved polymers in a chemically polydisperse sample depends on the chosen solvents via the interactions between the polymers and solvents.Using these distributions and average features of the dissolved and undissolved polymer, fractionation steps can be optimized.Furthermore, resulting distributions in chain length, functional groups, or other varying molecular features of polymers can be included in the same manner for subsequent steps in a production process.In the next section, we combine a distribution in polymer chain lengths and the number of functional groups per segment for a system consisting of lignin dissolved in various solvents.

The modeling of the solubility of lignin
In the previous section, we applied extensions of Flory-Huggins theory, including size and chemical dispersity, to hypothetical systems containing only variation in polymer size or the number of functional groups per segment.These computations can be applied to polymers that are physically and/or chemically polydisperse.Here, we apply these methods to the practical example of lignin dissolved in a range of solvents, using experimental measurements of the molecular mass distribution and number of functional groups per segment.
The theoretical calculations in this section use input from earlier work [66].These experiments used 94% pure Protobind 1000 lignin and a mixture of methanol and ethyl acetate.Molecular masses were measured using GPC and the concentrations of functional groups using NMR For more details we refer to [54].Based on the mass and functional groups distributions, we expect a log-normal distribution for the sizes of the lignin molecules and propose an exponentially decreasing volume fraction of functional groups with the chain length of the functional groups.This can be described by: Since θ k decreases exponentially with k and M k scales linearly with k, this also suggests that the number of functional groups per polymer segment decreases exponentially with the polymer chain length.A correlation between the number of functional groups per segment and the molecular size was earlier already observed by Sevastyanova et al [74].The dependency of θ k on k in equation (20) describes how functional groups are distributed over polymers with different lengths.The fitting procedure to find the fitting parameters is discussed in our earlier paper [66].The fitting parameters in equation ( 20) are e µ ≈ 2.7 • 10 3 , e σ = 1.2,  7).The parameter χ f represent here the interaction energy between each hydrophilic/hydrophobic connection, between polymer segments, solvent molecules and polymer segments and solvent molecules.Both forms of polydispersity are described with distributions, each with accompanying parameters.Most parameters are either from literature, related to specific measurements, or related to chemical features of lignin.Also, the number of properties described simultaneously provides extra confidence for the values of the parameters.At least one measurement of the molecular mass and function group distribution is necessary to make predictions for the solvent fractionation of a polymer sample.Though estimating or measuring these parameters requires effort, once these parameters are determined, they can be applied to a wide variety of situations in combination with readily available Hansen solubility parameters of solvents in the literature.A prediction of the solubility of a different type of lignin in a different solvent mixture is possible with only minimal changes in these parameters [66].
The theoretical average molecular mass and number of functional groups per segment of the dissolved lignin, can be determined using: Using the undissolved fractions, the average features of the residue can be determined.In the experiments, [66] NMR was used to measure the concentration of aliphatic, aromatic, carboxylic hydroxyl content.Only the total hydroxyl content was used for the theory.This neglects some of the complexity of the system, but reduces the number of fitting parameters.
In figure 5, the resulting theoretical predictions (curves) are compared with experimental data.Using these parameters, the trend in both the weight average molecular mass and the average number of functional groups per segment are close to the experiments.The minimum weight average molecular mass for the undissolved lignin and the maximum weight average molecular mass for the dissolved lignin match the total average molecular mass.At the same time, the maximum hydroxyl content of the undissolved lignin and the minimum hydroxyl content of the dissolved lignin are equal to the average total hydroxyl content.Only a small amount of lignin dissolves for large RED, minimally affecting the average values.When the RED is smaller than the RED where all lignin dissolves, the average dissolved molecular mass and hydroxyl content do not change anymore.The average molecular mass for the undissolved lignin approaches infinity at this value.
Also the experimentally observed and theoretically described yield of solubilized lignin can be compared.This yield can be calculated as a function of the RED using equation (15).In this equation, equation ( 11) is used as input for q.In equation (11), we use both limits and the averaging equation ( 7) and the equation for θ s from equation (7).Finally, we inserted the parameters in the resulting form of equation ( 15) and solved the equation.With η known, equation ( 13) can be used at values for k while varying the RED to find the  dissolution curves of individual chain lengths of lignin, or varying k while keeping the RED constant to determine the theoretical molecular mass distributions.
Figure 6 shows the individual amounts of lignin dissolved and the total fraction of lignin that dissolves as a function of the RED.For individual polymer chain lengths, the amount of polymer dissolving increases from close to zero to its maximal value in a small range of RED-values.The combination of all these curves results in a much wider range of RED-values, that corresponds better to the experimentally measured dependency.The agreement between theory and experiment for the solubility for a different type of lignin in a different solvent mixture provides extra support for this correspondence [66].The addition of the predictions for the average molecular features presented here shows that using the theory described in the previous sections simultaneously predicts the RED-dependency of the yield, weight average molecular mass, and hydroxyl content of the dissolved lignin.
Figure 7 shows the theoretical and experimental molecular mass distributions of the dissolved part of lignin.The intensities are normalized using the maximal measured intensity of the measurement with the highest yield.Although the shapes are not exactly equal, there are trends in the shape that can be observed in both the theoretical and experimental results.In both the theoretical and experimental results, the maxima of the distribution shift to smaller dissolved fractions.Also, both results suggest that mainly the higher molecular mass polymers stop dissolving at higher RED.Together with the results for the yield, average molecular mass, and hydroxyl content, this supports that the theory can be used to include size and chemical variation distributions in Flory-Huggins theory calculations to describe multiple characteristics of dissolving polymers simultaneously.

Conclusions
In this paper, it is shown that adaptations of Flory-Huggins theory can be used to describe solvent fractioning of polymers.A complication in using bio-based or recycled feedstock for polymeric materials is the variation in material properties due to the variation in the molecular features of the involved polymers.Comparing monodisperse and polydisperse calculations shows that the degree of polydispersity in chain length and the variation of functional groups affect the characteristics of polymer samples.Solvent fractionation improves the reproducibility of these properties by selecting a subset of molecular features.It is demonstrated here how the molecular features of solubilized lignin can be described using polydisperse Flory-Huggins theory relations.The presented relations make it possible to predict how chemical and physical dispersity affect polymer fractionation in different solvents.Comparisons with experiments shows simultaneous quantitative agreement for the yield, average molecular mass, and hydroxyl content and qualitative agreement for the molecular mass distributions, by including polydispersity in the calculations.

Figure 2 .
Figure 2. Volume fraction (a) and yield (b) of a monodisperse polymer when 3 vol% of polymer is added to a solvent.The different dashed curves represent polymers with different chain lengths.The colored regions represent the range of possible valuesfollowing from the approximations in the maximal and minimal limit as discussed in section 3.1.The curve for Mp = ∞ is the dependency obtained directly from the Hansen solubility parameters (η = 1 for RED < 1 and η = 0 for RED > 1), which is valid in the limit of infinitely long polymers.The gray curves represent the numerical results of the Flory-Huggins theory equations (4) and (5).

Figure 3 .
Figure 3. Yield as a function of the RED, for log-normal distributed chain length polydisperse polymers (a) and polymers with a normal distribution in the number of functional groups per segment (b).The average polymer chain length (⟨M⟩ = 100) and interaction energy (⟨θ⟩ = 1 3 ) are kept constant, while the width of the (log-)normal distribution σ is varied.The top and bottom curves give the solutions in the theoretical limits and Vp Vtot = 3 100 .

Figure 4 .
Figure 4. Distributions of the dissolved and undissolved fraction of polymer at four RED values.Parameters correspond to the parameters used for figure 3, with for the log-normal chain length distributions σ = 1.00(a) and for the normal interaction energy distributions σ/µ = 0.30 (b).

Figure 5 .
Figure 5.The weight average molecular mass and weight average number of functional groups per segment from equations (20) and (21) (curves) and experimental data[66].The solid curves show the features of the dissolved lignin and the dashed curves of the undissolved lignin.The error bars for the experimental data in the presented figures are standard deviations, with each experiment performed twice.

Figure 6 .
Figure 6.Theoretical individual and total dissolved amounts and experimental dissolved amounts of lignin as a function of the RED.The individual curves represent the amounts dissolved of individual species k, with k increasing with 50% for every next curve (rounded up, starting with k = 1).The area between the theoretical curves gives the range of possible values following from the minimal and maximal assumptions in equation (11).

Figure 7 .
Figure 7. Theoretical and experimental molecular mass distributions of dissolved lignin in varying solvents.The experimental curves correspond with the experimental data shown in figures 5 and 6.The equally colored curves in both figures correspond to the same RED-values.