Effects of Polydispersity on the Phase Behavior of Additive Hard Spheres in Solution

The ability to separate enzymes, nucleic acids, cells, and viruses is an important asset in life sciences. This can be realised by using their spontaneous asymmetric partitioning over two macromolecular aqueous phases in equilibrium with one another. Such phases can already form while mixing two different types of macromolecules in water. We investigate the effect of polydispersity of the macromolecules on the two-phase formation. We study theoretically the phase behavior of a model polydisperse system: an asymmetric binary mixture of hard spheres, of which the smaller component is monodisperse and the larger component is polydisperse. The interactions are modelled in terms of the second virial coefficient and are assumed to be additive hard sphere interactions. The polydisperse component is subdivided into sub-components and has an average size ten times the size of the monodisperse component. We calculate the theoretical liquid–liquid phase separation boundary (the binodal), the critical point, and the spinodal. We vary the distribution of the polydisperse component in terms of skewness, modality, polydispersity, and number of sub-components. We compare the phase behavior of the polydisperse mixtures with their concomittant monodisperse mixtures. We find that the largest species in the larger (polydisperse) component causes the largest shift in the position of the phase boundary, critical point, and spinodal compared to the binary monodisperse binary mixtures. The polydisperse component also shows fractionation. The smaller species of the polydisperse component favor the phase enriched in the smaller component. This phase also has a higher-volume fraction compared to the monodisperse mixture.


Introduction
The ability to separate enzymes, nucleic acids, cells, and viruses is an important asset which can be realised by using spontaneous asymmetric partitioning over two aqueous macroscopic phases [1]. The formation of such two phases can be induced by mixing two different types of macromolecules in an aqueous phase and depends on the polydispersity of the macromolecules. The prediction of such phase formation thus forms an important aspect of life science applications. Apart from applications, the formation of separate aqueous phases from the cytoplasm has also received interest in the last decade [2]. Interestingly, the pre-assembly mechanism during evolution as described recently [3] may be speculatively related to the same separation processes.
The formation of separate aqueous phases is usually studied using mixtures of two monodisperse components. However, components in nature are usually not that simple or well-defined. Often, components will exhibit polydispersity in terms of their size, shape, and charge, which is often ignored when studying phase behavior. In order to study the effects of polydispersity theoretically, we chose the following model system: hard spheres. Their separation into two phases is driven by two different physical mechanisms. One mechanism involves only excluded volume interactions, where the minimum distance between the particles is determined by the sum of their respective radii [4], where only a certain asymmetry in the sizes of the particles in the mixture is necessary [5]. This asymmetry leads to the depletion of small spheres around the large spheres and as a result to an effective attraction (depletion interaction) between the larger spheres [6]. This mechanism is referred to as additive hard sphere (HS) model. The case resembling the second physical mechanism is characterised by the minimum distance between the particles being larger or smaller than the sum of their respective radii. This case is referred to as nonadditive hard sphere (NAHS) model. In both cases, upon phase separation, the mixture will demix into two (or more) phases, each enriched in one of the components. In this work, we will focus on the first type: binary mixtures with significant asymmetry in their size.
In most studies on the phase behavior of binary mixtures, the polydispersity of the components is ignored. However, from experiments with for example gelatin and dextran, it is found that polydispersity has an influence on the phase behavior. The polydispersity of both components leads to significant fractionation, especially for dextran [7].
Polydispersity has an effect on the depletion interaction. With increasing polydispersity, the repulsive barrier decreases, leading to an enhanced rate of flocculation of the large colloidal particles [8]. Additionally, Waltz and co-authors [9] studied the depletion interaction in a solution of normally distributed macromolecules and showed that polydispersity can lead to increased flocculation through the formation of secondary potential energy minima. In later studies, they found that polydispersity significantly lowered the magnitude of the repulsive structural barrier, which can be understood in terms of a change in the depletion of the macromolecules from the gap [10]. The polydispersity of the smaller component affects the pair potential between the large particles [11]. This effect on the depletion interaction can stabilize the particle suspension in the short term but will still destabilize over time [12].
Studying the phase behavior of polydisperse mixtures is challenging, since a polydisperse component effectively consists of a large number of sub-components, each with a different size and possibly also different shape or charge [13]. Some theoretical work has been conducted on predicting the phase behavior of polydisperse components. Cotterman and co-authors used continuous and semi-continuous distributions to predict the fluid-vapor phase diagram of polydisperse components [14,15]. Santos and co-authors [16] studied the phase behavior of polydisperse compounds such as polystyrene and polyethylene glycol. They found that the polymer polydispersity played a crucial role in the phase behavior: the broad size distributions lead to a wide range of depletion attractions, giving rise to spinodal decomposition and preventing gelation. Bellier-Castella and co-authors [17] used a van der Waals approximation for free energy to study the phase behavior of polydisperse fluids composed of spherical particles. They found the onset of a three-phase co-existence at a higher polydispersity. Fasolo and co-authors [18] studied theoretically the equilibrium phase behavior of mixtures of polydisperse hard-sphere collids and monodisperse polymers based on the Asakura-Oosawa model. They found that with polydispersity, significant fractionation occurred. Polydispersity delayed the onset of both gas-liquid and fluid-solid separation. Additionally, Sear and co-authors [19] used the Asakura-Oosawa model to predict the phase behavior between a monodisperse colloid and a polydisperse polymer. They found that polydispersity increased the extent of the fluid-fluid co-existence. Warren [20] studied the interaction between hard spheres with a bimodal size distribution and found that demixing caused additional size partitioning and fractionation. We note that an approximation up to the third moment distribution used in that work is equivalent to the second-order virial approximation. See also, for example, [21]. Kang and co-authors used a universal quasi-chemical (UNIQUAC) model to predict the phase behavior of aqueous polymer systems. They found that the polydispersity of the polymers enlarged the two-phase region considerably near the plait point and resulted in smaller miscibility regions far from the plait point. They also found that the average molecular weights of polymers in the phases differed significantly and these differences increased with a larger polydispersity; due to this fractionation, the polydispersity of each polymer was smaller in each child phase compared to the parent mixture [22]. Others modelled the phase behavior of non-additive hard-sphere systems using Monte Carlo simulations [23]. They found that, with increasing polydispersity, the miscibility region decreased and that the critical point shifted towards lower pressures. Additionally, Stapleton and co-authors used Monte Carlo simulations to predict the phase behavior of mixtures with fixed or variable polydispersity [24]. They found that mixtures even with a very small degree of polydispersity resulted in differences in the phase separation and the fractionation between the coexisting phases.
In this study, we aim to gain a better understanding of how size polydispersity influences the liquid-liquid phase behavior in binary mixtures, mainly on the position of the phase boundary, the spinodal, and the critical point. Next, we aim to predict the fractionation of the polydisperse component between the phases. We model the interactions between the different components using the second virial coefficient. We note that the validity of the second-order virial coefficient approximation is limited to dilute systems for hard spheres, but depending on the type of system, as in the aqueous polymeric systems in which we are interested for the above-mentioned practical applications, the validity may increase to the polymer overlap concentration and beyond [25]. This warrants the exploration of the second-order virial coefficient approach. The resulting calculations form a stepping stone for expansion to non-additive interactions, polydispersity effects, and multi-component mixtures for practically relevant systems.
We start the theoretical considerations by reviewing the interaction in a simple system of a solute in a solvent (Section 2.1). In Sections 2.2 and 2.3, we expand the second virial coefficients for solutions with one type of solute component to solutions with multiple distinguishable types of solute components. Section 2.4 describes the theory about the stability of a mixture, Section 2.5 describes the theory about the critical point, and finally Section 2.6 describes the theory about the phase boundary. We chose to first describe the existing theory in order to more easily explain the expressions we used in our calculations. With the expressions in Section 2, we calculated the phase behavior for different mixtures with varying polydispersity. In Section 3, we discuss the resulting phase diagrams, in Section 3.1 we divide the polydisperse component into two sub-components, and in Section 3.2 we increase the number of sub-components to nine. Finally, we discuss the fractionation of the polydisperse component in Section 3.3.

Theory
We start by deriving the equations of state for dilute solutions. Next, we derive the virial expansion for solutions with one solute component. Subsequentl, y we derive the second virial coefficient for mixtures with an arbitrary number of distinguishable components. This gives us all the parameters we need to define the stability boundary, the critical point, and phase boundary of a mixture. The resulting system of equations is solved in Matlab R2017b.

Dilute Liquid Solutions
We consider a two-component solution, in which one component is the solvent and the other component is the solute. We define N s as the number of solvent particles and N ν as the number of solute particles in a volume V at a temperature T. The total number of particles in the system is then N = N ν + N s , and since we assume a dilute solution N s >> N ν . The system is in constant thermal contact with the environment, and both the volume and the number of solute and solvent particles are fixed (canonical ensemble) [26].
The sum of the kinetic (K) and potential energies (U) of the system represents the Hamiltonian (H) of the system, given by: in which: where p i is the impulse of particle i, m ν is the mass of a solute particle, m s is the mass of a solvent particle, φ ij is the pair potential between particle i and j, and r i is the position of particle i. The canonical partition function (Z) describes the statistical properties of the system for a given temperature, volume, and number of particles. The partition function is the sum of all the different individual energy states in which the system can exist. The states of the system are specified by both the position and the momentum of the particles. Applying the partition function to dilute solutions makes it possible to reduce the many-body problem in statistical mechanics to problems of one-body, two-body, three-body, etc.
where h is Plank's constant and β = 1 kT , in which k is Boltzmann's constant.
Other thermodynamic variables, such as the Helmholtz free energy, the pressure and the chemical potential can be expressed in terms of this function or its derivatives. The Helmholtz free energy (A) for this system is then given by [26]: With the differential of the free energy given by: Therefore, the pressure (p) is given by: and the chemical potential (µ i ) for component i is given by: Since we focus on particles with hard sphere interaction, we can integrate out the momentum integrals in Equation (4).
where Λ = h 2 2πmkT 1/2 is the mean thermal wavelength and Q the configuration integral.
The configuration integral is the integral over all possible configurations of the N molecules in the system: The first three configuration integrals are: The configuration integral Q 1 indicates that there is only one particle present in our selected volume. Q 2 indicates there are two particles present in our selected volume: interacting or not interaction. Q 3 indicates there are three particles present in the volume. These particles can interact with each other or not interact. The number of combinations of interactions increases significantly when increasing the number of particles.
The configuration integrals can be represented in diagrams where each dot is a particle present in the volume and a line between the dots indicates interaction: In our analysis, we consider solutions where the solvent particles are present in a much larger number than the solute particles. This means that solute particles have relatively low influence on the statistics of the solvent particles. Following the MacMillan-Mayer theory we can describe the interactions between the solute particles by a potential of mean force equation and thus can apply additivity of the particle-particle interactions in Equation (3) [27].
The potential of mean force for dissolved particles (W) is defined according to: Using Equations (11) and (19), the configuration integral becomes [28]: The Helmholtz free energy Equation (5) of the system then becomes the sum of the Helmholtz free energy of the solvent and the Helmholtz free energy of the solute:

Second Virial Coefficient of a Dilute Solution with a Single Solute Component
Similar to the expansion of the universal gas law by a virial expansion for real gasses, we can write a virial expansion for the osmotic pressure, Π, of a solution according to: With ρ as the number density of the component N ν V , B 2 as the second virial coefficient, and B 3 as the third virial coefficient. The second virial coefficient accounts for the increase in osmotic pressure due to particle pairwise interaction. The third virial coefficient accounts for the interaction between three particles. The equation can be expanded for higher densities with B n , the nth virial coefficient, which accounts for the interaction between n different particles. Until now, we have been using the canonical ensemble to describe the system. In the canonical ensemble, the number of particles (N ν + N s ), the temperature (T), and the volume (V) are fixed. The restraint of constant number of particles becomes tedious when accounting for the interaction between the different particles, therefore we will use the grand canonical ensemble to derive the virial coefficients. In the grand canonical ensemble, the temperature (T) and the volume (V) are fixed, as well as the chemical potentials (µ i ).
We can write the grand canonical partition function independent on N ν and N s by performing a transformation of the number of compounds by the chemical potential µ ν and µ s [26].
The equation of state for the system is given by: In which p is the pressure in the solvent reservoir and the osmotic pressure is given by: We can now define the activity as z ≡ e βµ Λ 3 The osmotic pressure can be written in terms of the logarithm of the grand canonical ensemble.
This can be written as: The coefficient b l is also known as the cluster-integral and it indicats interaction among l compounds.
Just as that configuration integrals can be written in diagrams (see Equations (15)- (17)), the cluster integrals can also be written as diagrams [29]: Two particle interactions can be defined by the Mayer f-function: This can be represented by: Substituting the Mayer f-function into the cluster integrals we get: In order to find the relationship between cluster integrals and the virial coeffcients, we need to do several substitutions and inversions using also (28): Inverting the series we obtain for the activity [30]: This can be substituted in the equation for the pressure (Equation (28)): Comparing this equation to Equation (22), we see that the second virial coefficient is equivalent to −b 2 : For hard spheres, we define the interaction potentials as for components with a diameter σ: The second virial coefficient then becomes:

Second Virial Coefficient of Solutions with Multiple Solute Compounds
When there are more distinct compounds (ν 1 , .., ν n ) in the solution with n the total number of distinguishable compounds, or species, there are two main types of two particle interactions that can occur: • interactions between indistinguishable components-i.e., components of the same species: • interactions between distinguishable components, i.e., components of different species: We can write the configuration integral Q N in general as: In which: N = n ∑ i N ν i or the total number of particles in the configuration and x and y can be of any type ν n in the mixture. The general equation for the partition function in the grand canonical ensemble then becomes: In the case of two-particle interaction we have interaction between components of the same species and interaction between components of different species, so we obtain for Q 2 two types of configuration integrals (comparable to Equation (13)): In which x and y can be of any type ν n in the mixture and y = x. Additionally, for the cluster integrals we obtain two types of integrals (comparable to Equation (36)): In which x and y can be of any type ν n in the mixture and y = x. The general equation for the osmotic pressure then becomes (comparable to Equation (38)): For the second virial coefficient, we obtain two types (comparable to Equation (41)): In which x and y can be of any type ν n in the mixture and y = x. Note: we define a B * to have all the second virial equations of the same form: −2π ∞ 0 drr 2 f (r), with f (r) dependent on the type of interaction.
For additive hard sphere interaction, the interaction potential for particles of different species is given by: With σ ij = (σ i + σ j )/2, the distance between the centers of the two components. When the interaction is not additive, the distance of closest approach of the centers of the two components becomes: σ ij = 1/2(σ i + σ j )(1 + ∆), in which ∆ accounts for the non-additivity of the interaction between the particles that are different.
The general equation for the osmotic pressure for a dilute mixture is then given by: The second virial coefficients can be represented in matrix form: In which we abbreviate the notation B ν i ν j to B ij , and similarly, the densities ρ ν i by ρ i .

Stability of a Mixture
The stability of a mixture is dependent on the second derivative of the free energy. If the second derivative of the mixture becomes zero, the mixture is at the boundary of becoming unstable. Unstable mixtures have a negative second derivative [31,32].
The free energy of a mixture is given by [26]: The differential is given by: In which the chemical potential (the first partial derivative of the free energy with respect to number of particles (N i )) for component i is given by: For a mixture with n distinguishable components, this second partial derivatives can be represented by a n × n matrix of the first partial derivatives of the chemical potential of each component.
This results in the following general stability matrix: When this matrix is positive definite, the mixture is stable [33]. Based on this criterion, when one of the eigenvalues is not positive, the mixture becomes unstable. When the matrix has one zero eigenvalue and is otherwise positive definite, the mixture is on the spinodal and is at the limit of stability [34].
In the case of a binary mixture (n = 2), the spinodal is also equal to the determinant of matrix M 1 . For mixtures with more components, this is not always the case anymore, as with an increasing number of components, there are an increasing number of eigenvalues for matrix M 1 that can become zero [33]. This can be resolved by checking if the stability matrix is positive definite for small changes of in the concentrations of the components near the concentration where the det(M 1 ) is zero.
The determinant of matrix M 1 for a monodisperse binary mixture is given by: Often, however, components are not a 100% monodisperse. Let's now investigate how the equation for the spinodal of the mixture changes when we introduce polydispersity in one of the components. We define a binary mixture in which one component (component 2) is polydisperse. The concentration of each of the particles in the polydisperse component can be represented by: With: Then: Each of the components in this mixture has a corresponding virial coefficient and cross virial coefficient.
Let us investigate the equations for a simple polydisperse mixture, in which the polydisperse component consists out of two sub-components (a and b, n = 3) (Figure 1). For the density and for the virial coefficient matrix we obtain: The stability matrix becomes: This results in the following determinant for the stability matrix: By increasing the number of sub-components, the number of terms in this determinant increases rapidly. This forms an incentive to try and treat the polydisperse component as if it is effectively one component. A natural and convenient choice for this route is coupled to the experimental determination of virial coefficients using membrane osmometry [35]. Namely, membrane osmometry yields values that are number averaged. Thus, we choose number-averaged virial coefficients.
The number averaged virial coefficient of a mixture can be written as: In which B ii is the second virial coefficient of the ith particle, B ij is the second cross virial coefficient of the ith particle and the jth particle, and x i is the fraction of the ith particle, ∑ x i = 1.
Using this definition, we can map the polydisperse mixture by a 2 × 2 matrix of virial coefficients. We will refer to this 2 × 2 matrix of effective virial coefficients.
B 11 e f f = B 11 For the mixture we considered in Equation (64), the effective virial coefficients become: The effective stability matrix for this mixture then becomes: 2B 11 e f f + 1 The determinant then becomes: It is clear that Equations (66) and (71) are different. Using the effective virial coefficients to determine stability of the mixture possibly results in deviations.

Critical Points
In a binary mixture, the critical point is a stable point which lies on the stability limit (spinodal) [34] and where the phase boundary and spinodal coincide. In mixtures of more components these become plait points. Critical points and plait points are in general concentrations at which two phases in equilibrium become indistinguishable [36].
There are two criteria that have to be used to find critical points. The first one is det(M 1 ) = 0, which is the equation for the spinodal. The other criterion is based on the fact that at the critical point, the third derivative of the free energy should also be zero. For a multicomponent system, this criterion can be reformulated using Legendre transforms as det(M 2 ) = 0 [31,37], where: Matrix M 2 is matrix M 1 with one of the rows replaced by the partial derivatives of the determinant of matrix M 1 . Note: it does not matter which row of the matrix is replaced.
For a monodisperse binary mixture, this results in the following two matrices for the critical point: In addition: The set of equations that needs be solved for the critical point is: For the earlier considered polydisperse mixture containing the two sub-components (a and b, n = 3), we obtain: With: In general, P i can be found using the following equation: In which M 1,(ii) is the minor of matrix M 1 at the i th -row and i th -column.
Combining det(M 1 ) and det(M 2 ) results in the following set of equations: Since, in this set of equations, there are more higher-order terms present, it is possible that this results in multiple plait points, depending on the concentration of each of the components in the mixture. Care should be taken that the solutions of the set of equations have concentrations at the limit of stability; this can be done by checking the eigenvalues of the stability matrix.
If we use the effective virial coefficients for this mixture as defined in Equation (69), we obtain: The determinant of this matrix becomes: This results in the following system of equations for the critical point: Additionally, for the third derivative of the Helmholtz free energy, we see that reducing the polydispersity by using the effective virial coefficients, results in fewer terms in the equation and possible deviations in determining the critical point.

Phase Boundary
When a mixture becomes unstable and phase separates into two or more phases, the chemical potential of each component and the osmotic pressure is the same in all phases [26].
where the phases are denoted by I, I I, .... For a system that separates into two phases, we obtain, using Equations (56) and (60), the following set of equations for the general case of n components: This set of equations has 2 × n unknowns and n + 1 equations. The set of equations can be solved by fixing one of the concentrations for one phase and the ratio of the concentrations of the other components for the same phase. To solve this set of equations in order to find the concentration of each component in each phase, without fixing any of the concentrations, we need therefore an extra set of equations.
This extra set of equations stems from the fact that during phase separation no particles are lost and no new particles are created. The total number of components in the system is therefore given by: Additionally, the total volume, V, of the system does not change. The total volume of the system is given by: The concentrations of each component in each phase are thus given by: The total number of compounds in the system can be found using: This can be rewritten to: This results in an extra set of n equations and one more unknown (α). The complete set of equations to solve for the binodal then becomes: . .
Systems for which there are more than two distinguishable components (n > 2) can theoretically have more than two coexisting phases, according to the Gibbs phase rule. With an increasing number of phases, the set of equations to solve increases as well. The number of equations needed to solve for an arbitrary number of f phases is: f × n + f − 1. The set of equations then becomes: . . With: Now that we have the general equations for the phase boundary, let us investigate the set of equations we need to solve for the mixtures we defined earlier. For a monodisperse binary mixture the set of equations to solve for the phase boundary is given by: ln(ρ I 1 ) + 2B 11 ρ I 1 + 2B 12 ρ I 2 = ln(ρ I I 1 ) + 2B 11 ρ I I 1 + 2B 12 ρ I I 2 ln(ρ I 2 ) + 2B 12 ρ I 1 + 2B 22 ρ I 2 = ln(ρ I I 2 ) + 2B 12 ρ I I 1 + 2B 22 ρ I I For the polydisperse binary mixture we considered earlier (with two sub-components a and b, n = 3) this set of equations becomes: Note that the ratio between ρ I 2 a and ρ I 2 b is not necessarily the same as the ratio between ρ I I 2 a and ρ I I 2 b , since fractionation between the components can occur. Using the effective virial coefficients, we obtain: It is clear from the equations that when using the effective virial coefficients for calculating the phase boundary, all information about the (changes in) distribution of the polydisperse component becomes untraceable.
We are calculating the transition from one-phase to two-phase systems and approach this point from the one-phase lower-concentration regime.

Results and Discussion
In this work, we calculated the liquid-liquid phase diagram for a variety of binary additive mixtures of a small hard sphere A and a larger hard sphere B with a size ratio q = σ A /σ B = 1/10. We started by calculating the phase diagram of this monodisperse mixture ( Figure 2) and gradually introduced polydispersity into the composition of component B (Figures 3-6). Component B is characterized by a degree of polydispersity (PD), defined by: For all particles, the concentrations are expressed as a dimensionless parameter according to η = πρσ 3 6 . We calculated the critical point, the phase separation boundary, and the spinodal of the various mixtures. Next to that, we also investigated the composition of the child phases, volume ratio between the phases (α), and the fractionation of the polydisperse component B for a specific parent mixture. In order to be able to compare our results with the monodisperse case, we first refer to Figure 2, where we see that the binary mixture phase separates at very low volume fractions of component A (η A crit = 0.007) and significantly higher concentrations of component B (η B crit = 0.267). This high asymmetry in the position of the critical point and phase boundary has also previously been reported by [38]. Note that even though the volume fraction of the smaller spheres are very low, their number concentration is significantly higher than the large spheres.

Polydisperse Mixtures with 2 Sub-Components
In Figure 3, we show the phase diagram for the case of a slight polydispersity in component B. Component B consists of two sub-components and has a PD = 4.00. These components are additive spheres in two sizes (both present in the same amount), with the number average size of the mixture equal to that of the size of the monodisperse mixture of Figure 2. The mixture therefore consists of three components. We calculated the phase diagram using both the simplified 2 × 2 effective virial coefficient matrix described in the theory (we refer to this as the effective mixture) and the full 3 × 3 virial coefficient matrix (to which we refer as the polydisperse mixture). The difference between the phase boundary, spinodal and critical point of the monodisperse mixture and the effective mixture is negligible. We see however that the introduction of the polydispersity causes the critical point to shift to a higher volume fraction of component B (η A crit = 0.007, η B crit = 0.280) and that especially at lower volume fraction of component B the phase separation boundary shifts towards slightly lower packing fractions.  For Figure 4, we increased the size difference between the smaller and the larger spheres for component B, with the standard deviation twice the standard deviation of the spheres in Figure 3 and PD = 8.00. The patterns we saw in Figure 3 are more pronounced for this mixture: the increase in size difference causes the critical point to shift to higher packing fractions for component B and lower packing fractions for component A (η A crit = 0.006, η B crit = 0.315). The phase boundary, i.e., the binodal, shifts to significantly lower packing fractions, especially at lower volume concentrations of component B. For this mixture we see also that the spinodal of the polydisperse mixture shifts towards lower packing fractions of A. There is a slight difference between the positions of the binodal, spinodal, and critical point of the monodisperse mixture and effective mixture.
In Figure 5, we introduced skewness in the size distribution of component B. For the first two mixtures (Figure 5a,b), the ratio between the bigger and the smaller sub-component was 25/75. The polydispersity for both mixtures is the same, PD = 6.93. For the other two mixtures (Figure 5c,d

Polydispersity with 9 Sub-Components
Now that we have a bit of an understanding of how polydispersity influences the critical point and the phase boundary, we increase the number of sub-components of B. The mixtures in Figure 6 thus consist of ten components (one component A and nine components B in varying amounts and sizes, with different degrees of polydispersity, dependent on the considered distribution). We again calculated the phase diagram using both the simplified effective 2 × 2 virial coefficient matrix (the effective mixture) and the full 10 × 10 virial coefficient matrix (the polydisperse mixture). Table 1 gives the critical points for the polydisperse mixtures and Table 2 allows for an easy comparison of the different distributions.
From the figure and table, we can conclude that the standard deviation of the polydisperse component plays a big role in moving the critical point and phase boundary. The standard deviation for B in Figure 6b is twice the standard deviation for B in Figure 6a. It is also clear that the type of distribution plays a significant role in the concentration of the critical point and position of the phase boundary. The sizes in B for Figure 6b,e are the same, however, each size is present with a different frequency. The distribution in Figure 6b is Gaussian and the sizes in Figure 6e are bimodal, which means that the particles with sizes just larger and just smaller than the mean are present in a larger number. This causes to shift the critical point to slightly lower packing fraction of B in Figure 6e compared to Figure 6b. The distribution of B in mixture Figure 3 is comparable to the mixture in Figure 6a, just with fewer sub-components. The position of the critical point is for both mixtures very comparable. In the same way is the distribution of B in Figure 4 comparable to the mixture in Figure 6b,e. However, the increased number of sub-components results in a slightly higher and lower concentration of the B component at the critical point respectively.
The distributions of B in mixtures in Figure 6c,d are skewed. In this regard they are comparable to the mixtures in Figure 5. For these mixtures we see that the right skewed distribution also causes the critical point to move to higher concentrations of B compared to the mixture with the left skewed distribution.

Fractionation of Polydisperse Component
Upon phase separation, particles will move to a preferential phase in order to minimize the Helmholtz free energy. One phase is enriched in component A, whilst the other is enriched in component B. Even though each phase is enriched in one component, the other component is still present in lower concentrations. We investigated the phase separation of a specific parent mixture (η A parent = 0.010, η B parent = 0.200) for the different distributions of component B, in terms of volume fraction of both components in each phase, degree of polydispersity of component B, average size of component B in child phase compared to the average size of component B in the parent phase and the volume fraction of the phases (α), see Table 2.     Figure 5a, we observe that the smaller sub-components favor the top phase (the phase enriched in the small particles A), while the larger sub-components favor the bottom phase (the phase enriched in the larger particles B). We like to note that these observations as obtained from the virial approach are in line with previous theoretical work using the UNIQUAC model by Kang and Sandler [22] and the Florry Huggins theory [39] and also experimental work [7,40,41]. Ref. [42] also note that the skewness of the polydisperse parent distribution plays an important role in the fractionation. Additionally, ref. [23] found that the largest polydipserse particle favored the phase poor in the smallest particle. This follows from the relatively increased size incompatibility between the smallest particle and the largest particle in the system.
Depending on the type of distribution, phase separation causes the degree of polydispersity to decrease. This is in both phases for all mixtures with a symmetric distribution of B in the parent phase. For these mixtures, the degree of polydispersity is at its lowest for the phase enriched in A.

Conclusions
We find that the largest species