Oligomer/Polymer Blend Phase Diagram and Surface Concentration Profiles for Squalane/Polybutadiene: Experimental Measurements and Predictions from SAFT-γ Mie and Molecular Dynamics Simulations

The compatibility and surface behavior of squalane–polybutadiene mixtures are studied by experimental cloud point and neutron reflectivity measurements, statistical associating fluid theory (SAFT), and molecular dynamics (MD) simulations. A SAFT-γ Mie model is shown to be successful in capturing the cloud point curves of squalane–polybutadiene and squalane–cis-polybutadiene binary mixtures, and the same SAFT-γ Mie model is used to develop a thermodynamically consistent top-down coarse-grained force field to describe squalane–polybutadiene. Coarse-grained molecular dynamics simulations are performed to study surface behavior for different concentrations of squalane, with the system exhibiting surface enrichment and a wetting transition. Simulated surface profiles are compared with those obtained by fitting to neutron reflectivity data obtained from thin films composed of deuterated squalane (d-sq)–polybutadiene. The presented top-down parametrization methodology is a fast and thermodynamically reliable approach for predicting properties of oligomer–polymer mixtures, which can be challenging for either theory or MD simulations alone.


■ INTRODUCTION
Quantitative prediction of the effect of oligomers in mixtures with polymers is a challenging problem and is the subject of considerable active research. 1,2 In many industrially relevant polymers, oligomers are ubiquitous, appearing as additives, (e.g., tackifiers, plasticizers, compatibilizers, and reactive species) and as low molecular weight components of many polydisperse compounds. 3,4 Understanding their role is therefore of critical interest in the development of new applications. 5−7 In typical oligomer−polymer mixtures, microand macroscopic properties are affected by the molecular weight distribution, the mixture composition, the current temperature, and thermal history, in addition to chemical effects arising from (for example) differences in polarizability, polarity, and the ability to form hydrogen bonds.
Oligomers play a significant role in influencing polymer properties both in the bulk and at an interface. In bulk polymers, slow phase separation of oligomers can prove a major problem when designing stable polymer mixtures for specific applications. 8 However, oligomers also exhibit a significant (sometimes beneficial) influence on polymer surfaces through surface migration, wetting transitions, 9 and surface segregation (including blooming). 10−15 Much progress has been made in studying phenomena such as oligomer phase separation, including theoretical approaches such as lattice integral equation theories. 16 Moreover, the role of entropy in causing surface migration of low molecular weight species is well-understood. 12,17 However, the influence of the surface on phase separation and the influence of bulk phase separation on surface migration, the formation of wetting layers, and surface blooming of additives are more complicated to predict quantitatively. 17 Coarse-grained (CG) simulation models 18−25 provide a way of studying the influence of oligomers at a molecular level. Here, accessible times scales can be many orders of magnitude longer than those available from atomistic simulations. CG models provide a way of studying surface enrichment 26 and provide a direct link to surface sensitive experimental measurements, such as neutron reflectivity. 2,27−30 However, there is a considerable challenge associated with the production of quantitative coarse-grained models that work well in both the bulk and at an interface and can properly represent the mixing thermodynamics (that drives phase separation) over a range of concentrations. For example, bottom-up approaches to develop coarse-grained models from an underlying atomistic description, such as iterative Boltzmann inversion (IBI) 31,32 and multiscale coarse graining (MSCG) 33 (force matching), have been extensively used for polymers. 22,23,26,34 However, these approaches usually rely on atomistic structural information from intra-and intermolecular distribution functions, typically at a single state point. Consequently, they often perform poorly in predicting mixing free energies that determine phase separation and are often not sufficiently transferable to represent accurately the effects of temperature and pressure changes in systems with chemically distinct interaction sites. 24 In contrast to bottom-up models, top-down approaches, based on SAFT theories, 35−37 and MARTINI-like models 25,38−40 focus on capturing macroscopic thermodynamic experimental information in the development of CG polymer models. Hence, they might be expected to perform better under conditions where oligomer−polymer phase separation may be present. The most recent incarnation of the SAFT equation of state, SAFT-γ Mie, 36,37 shows improved accuracy over previous SAFT-VR versions by including higher-order perturbations. It employs a Mie potential, which provides more adaptive models and which demonstrates an improved link to molecular simulations. 35,41,42 SAFT-γ Mie has successfully been applied to describe thermodynamic properties of low molecular weight compounds and their mixtures. 43−53 Excellent agreement with experimental results has been demonstrated even for secondorder properties such as heat capacities. 36,37,44 Few SAFT-γ Mie models have been developed for higher molecular weight compounds. However, two recent successful examples involve the prediction of polyethylene densities 37 and the study of polystyrene melts and polystyrene solutions in alkanes. 35 These studies suggest that it may be possible to use SAFT-γbased approaches to develop successful theoretical descriptions of oligomer/polymer systems and use the link between SAFT-γ Mie and molecular simulation to study surface behavior.
In this paper we use a combined SAFT-γ Mie and coarsegrained MD approach to build a coarse-grained model to describe a mixture of squalane and polybutadiene. This system acts as a simple model for a hot melt adhesive (HMA) in which molecular migration of the oligomeric species can take place. The theoretical approach correctly captures densities and the binary phase diagram and explains the changing surface concentration profile of squalane at the polybutadiene surface as a function of concentration. We validate the results by experimentally measuring cloud points of squalane in polybutadiene and by measuring vertical concentration profiles of deuterated squalane in a thin polymer film by neutron reflectometry (NR).
Phase Diagrams. Polymer films of ≈1 mm thickness, containing various proportions of squalane, were cast from 10% (w/w) toluene solutions in Teflon dishes. The solvent was evaporated at room temperature and pressure. Evaporation was considered complete when no more weight loss could be measured with a precision scale (typically more than 4 days). The use of a vacuum oven was avoided to ensure that a minimum amount of squalane would be lost during sample preparation. The films were peeled from the Teflon dish and transferred onto a glass slide. The turbidity of the samples was assessed by the naked eye at temperatures between 0 and 60°C to determine the cloud point, i.e., phase boundary. The Peltier stage of an AR2000ex rheometer was used to accurately heat the sample at a rate of 2°C/min. The cloud point was measured at least twice on each sample to ensure consistency of the results.
Neutron Reflectometry (NR). Oligomer/polymer thin films (≈100 nm) were prepared by spin-coating ∼1 mL of a 2% (w/w) toluene solution onto silicon blocks. To obtain the desired thickness, a Laurell spin-coater model WS-650-23 was used at a rotation speed of 2500 rpm (typical spinning time 30 s). The neutron reflectivity measurements were performed at the Rutherford Appleton Laboratory (ISIS pulsed neutron and muon source) in Chilton nr. Didcot, UK. The experiments used a standard air/layer/silicon setup for neutron reflectivity. 55 The d-sq/PB films were analyzed by using the SURF reflectometer, while for the d-sq/PBcis mixtures, the INTER reflectometer was used. The composition profiles were obtained for these films by measuring the specular reflectivity as a function of the momentum transfer Q = 4π sin(θ/λ), with λ the neutron wavelength and θ the incident angle. Generally, measurements are performed from before the critical edge (Q ≈ 0.01 Å −1 for silicon/air) to the point at which the signal is indistinguishable from the background (Q ≈ 0.25 Å −1 ). Such measurements typically lasted 2.5 h but were somewhat more rapid on INTER due to the large flux and simultaneous Q range. The analysis software IGOR pro with the Motofit package 56 was used to fit the neutron reflectivity data. Scattering length densities for the silicon and the native oxide layer were 2.07 × 10 −6 and ≈3.45 × 10 −6 Å −1 , which are consistent with values reported elsewhere. 2,57−59 Although it is necessary to fit NR data with model concentration profiles, the generic form of these has been previously verified by ion beam analysis, 2 and key features such as the total film thickness and thickness of the excess layer are unambiguously related to the frequency of fringes in the R(Q) data.

■ COMPUTATIONAL METHODS
SAFT-γ Mie Models. The statistical association fluid theory (SAFT) version used in this work is the SAFT-γ Mie 36,37 theory. In the nonassociating form of the theory, the Helmholtz free energy can be expressed as where A ideal accounts for the ideal kinetic and translational energy contributions, A monomer for the change in free energy due to attractive and repulsive interactions of the beads, and A chain for the free energy contribution by forming chains. The molecular model in the SAFT-γ Mie equation of state (EoS) is a chain of tangentially connected beads interacting through Mie potentials, which can be expressed for the interaction between bead i and j as Ä with i k j j j j j j y where ϵ ij is the potential well depth, σ ij is the bead diameter, and λ r,ij and λ a,ij are the repulsive and attractive exponents. For λ r = 12 and λ a = 6, the potential takes the usual form of the Lennard-Jones potential.
The unlike interactions can be obtained from the like interactions with the usual arithmetic combining rule for contact distances, σ ij = 1 / 2 (σ ii + σ jj ), and the geometric combining rule for the well depth Here, the additional parameter, k ij , measures the extent to which the unlike interaction well depth deviates from the values given by the standard combining rules. Similarly, for the repulsive exponent the additional γ ij parameter quantifies the deviation from the combining rule introduced by Lafitte. 36 For the attractive exponent, the same combining rule is applied, but without introducing a correction factor Our in-house SAFT program, 60 based on the equations in the work of Lafitte et al. 36 and validated against the results in this paper, and the gSAFT package (Process Systems Enterprise Ltd.) 61 were used to develop the SAFT-γ Mie model. The SAFT-γ Mie models of the PB/sq and cisPB/sq mixtures were developed in a two-step process. First, the pure component Mie potential parameters of PB, cisPB, and sq were optimized to reproduce densities across a wide range of temperatures and pressures (pVT data) 62−64 of each compound. The four Mie potential parameters ϵ, σ, λ r , and λ a were optimized simultaneously to minimize the average absolute deviation from the reference set of densities where N exp is the number of experimental reference points. Here, the experimental density, ρ exp , was compared with the density predicted by the SAFT-γ Mie model, ρ SAFT , for each data point k. The number of beads per chain was set to give two beads per monomer for each compound (see Figure 1), which corresponds to a bead mass of 27 g mol −1 for PB and cisPB and 35 g mol −1 for sq. This mapping was chosen to reduce the number of beads compared to the number of atoms (coarse graining). A further decrease of the number of beads was not feasible as it either required large repulsive exponents or could not capture the compressibility and thermal expansion with the same parameter set. Second, the unlike interaction parameters between PB and sq and cisPB and sq were adjusted to reproduce each experimental cloud point curve.
We note that for many low molecular weight systems λ a is often fixed at 6 in the fitting of SAFT-γ Mie models, as the exponents of optimal models are not completely independent. 65 In some of our fitting work we also tried fixing λ a = 6, which has the added advantage of reducing the number of variables to fit. However, this typically led to λ r greater than 12 and higher than those obtained with λ a free to vary. Higher repulsive exponents limit the time step that can be used in coarse-grained MD studies and thus restrict one of the major advantages of using coarse-grained MD. Hence, throughout this work, Mie parameters have been fitted for both repulsive and attractive exponents λ r and λ a .
In principle, SAFT can be extended to surfaces by combination with a theoretical method capable of representing the new spatial dimension while using the SAFT equations to represent the bulk free energy. Such methods are typically based on square gradient theory (SGT) 66−68 and classical density functional theory (DFT). 69,70 Here, however, we use a molecular dynamics (MD) approach to study surfaces based directly on the SAFT-γ Mie parameters obtained in this study, allowing surface segregation to arise naturally from the intermolecular interactions present in the model.
Coarse-Grained Molecular Dynamics Simulations. SAFT-γ Mie models have been used previously to develop force fields for a number of low-molecular-weight pure compounds 44−47 and mixtures. 48−52 Recently, such studies have been extended to lyotropic 53 and chromonic 25 liquid crystals and polymer solutions. 35 The latter were developed by fitting SAFT-γ Mie theory for component parts, i.e., building blocks, and in a second step, some or all of the unlike interactions were adjusted in the molecular dynamics framework. 35 However, in the current work we attempt to produce polymer and oligomer MD models directly from SAFT-γ Mie fits to experimental data.
The CGMD force field used in this work was based on the SAFT-γ Mie model developed above. Mie potentials were used to describe nonbonding interactions, together with LINCS bond constraints with equilibrium bond lengths of r 0 = σ ij , to provide tangentially bonded spheres of diameter σ ij (see Figure 2). No explicit angle and dihedral potentials were applied. However, nonbonded interactions were calculated for all sites with the exception of bonded neighbors, which were excluded.
All molecular dynamics simulations were performed using GROMACS 4.6.7. 71,72 The equations of motion were solved using the leapfrog integrator with a time step of 5 fs (4 fs for higher   If not stated otherwise, simulations were performed at a temperature of 298 K using the Nose−Hoover thermostat. The polymer chains and initial configurations were built with the Assemble! program. 73 Here, we restrict the maximum size of the polymer chain to 120 beads (Figure 2), as the full length of a 2 bead per monomer chain, matching the experimental molecular weight, would be prohibitively expensive to equilibrate due to slow diffusion, high viscosity, and slow chain relaxation. For pure squalane we simulate 181 chains (2172 beads). For pure polymers we simulate chains of length 30, 40, 60, and 120 beads (4200 beads in each case). We note that the longest (120 bead) polymer chain simulated is 10× longer than the squalane chain and is therefore expected to be long enough to capture, at least semiquantitatively, the behavior of squalane−polybutadiene mixtures.
For mixtures, the following procedure was adopted by using appropriate bead numbers for different weight fractions of squalane. A randomly generated vapor mixture configuration of 48 PB and 92 squalane chains was prepared initially (∼20% (w/w) squalane). The system of cubic box size was compressed at 500 bar and 450 K for 1000 ps and pre-equilibrated at 450 K and 1 bar for 500 ps by using the Berendsen thermostat and Berendsen barostat in the constant-NpT ensemble. A full equilibration run was then performed of at least 1.5 μs duration with the Nose−Hoover thermostat (coupling constant, τ T = 1 ps) and the Parrinello−Rahman barostat (coupling constant, τ p = 5 ps). The equilibrated cubic box was stacked 3 × 1 × 1 to give a 3× larger system size of ∼20K beads and further equilibrated. Additionally, a larger system size of 4 × 2 × 2 boxes with over 100000 beads was studied to look at the influence of finite size effects. The larger systems showed similar surface profiles as the smaller box of miscible systems. Lastly, a vacuum layer was added, and the system re-equilibrated in the canonical (constant-NVT) ensemble. Phase-separating systems, or systems at conditions close to the phase boundary, required extremely long equilibration times when started from a fully mixed initial state due to kinetically trapped layers. In these cases, elongated simulation boxes helped facilitate the formation of vapor−liquid interfaces.
Density profiles along the box were obtained at a series of time steps by using the GROMACS gdensity tool. As there will always be (at least) two interfaces, the slab was split into two, and the surface profiles were averaged over both interfaces. This was done by an automated postprocessing step for each saved coordinate set. Simulation snapshots were rendered with the VMD visualization package. 74,75 Finally, it should be noted that in both the simulations and SAFT work only single components or binary mixtures are considered, and so the effects of polydispersity are not modeled explicitly; however, the Mie parameters used come from a fit over the range of molecular weights present in the real polymer. In terms of phase boundaries this should be a reasonable approximation for cases, such as here, where the polymer chains are much longer than oligomer. However, in other cases, where a polymer contained a significant fraction of lowmolecular-weight oligomer as part of the molecular weight distribution, we would also expect to see surface enrichment of this species due to entropic reasons. 12 ■ RESULTS AND DISCUSSION SAFT-γ Models for Squalane, PB, and CisPB. Pure component SAFT-γ models for squalane, PB, and cisPB were obtained by simultaneously fitting the ϵ ii , σ ii , λ r ii , and λ a ii Mie potential parameters to reproduce the experimental liquid densities across a wide range of temperatures and pressures (Figure 3). The fits obtained are excellent, with average absolute deviations obtained of 0.016% for sq, 0.036% for PB, and 0.11% for cisPB, across a wide temperature (60 K for squalane, 175 K for PB, and 198.1 K for cisPB) and pressure range (400 bar for squalane and 2000 bar for PB and cisPB). The deviations of the model predictions are below the experimental uncertainties (e.g., 0.4% for ref 63). For each system, SAFT-γ Mie leads to a single transferable set of Mie parameters that can be used to represent the entire range of experimental data (Table 2).
Although a wide range of experimental state points is covered in Figure 3, we note that only PVT data have been considered in the fitting process. For many low-molecularweight molecules, additional high-quality data are often available: for example, derivative properties, such as speed of sound measurements, can often be measured to high precision and are considered challenging to fit for thermodynamics models. 76 While such data are not easy to obtain for polymers (e.g., for high-viscosity reasons in the case of speed of sound measurements), including additional high-quality thermodynamic-derived data in the fitting process may further improve upon the transferability of the Mie parameters obtained in this work.
Molecular Dynamics Simulations of Pure Components. Figure 4 shows a comparison of simulated and experimental densities for pure squalane and pure polybutadiene. The simulation results are surprisingly good and testify to the effectiveness of SAFT-γ Mie as a procedure for generating coarse-grained simulation models without reference to an underlying atomistic description. For squalane, agreement is particularly good at higher temperatures (<0.01% error at 343.15 K) but even at the lowest experimental temperature in Figure 4 of 278.15 K, where we would expect some differences between simulation and experiment due to the lack of chain rigidity (no angle bending terms are included), the error is only 0.4%.
For polybutadiene, we expect density to increase with increasing chain length for the moderate chain lengths simulated here but for the density to become largely independent of chain length as the chain becomes sufficiently long (i.e., an indicator of true polymeric behavior). It is usual to represent the specific volume, v, by a hyperbolic equation 77 Macromolecules pubs.acs.org/Macromolecules Article where N is the number of backbone carbon atoms in a chain (2 per coarse-grained bead), v ∞ is the specific volume at infinite chain length, and v 0 is a proportionality constant describing the rate of decrease of v with increasing chain length. For PB chains with 120, 60, 40, and 30 coarse-grained beads, eq 8 yields a straight line plot (inset diagram in Figure 4) allowing v ∞ to be estimated as 1.2202 cm 3 g −1 . This is an error of 2.1% in comparison to experiment for chains of M w = 280000 g mol −1 , ignoring any effects arising from polydispersity or the presence of 1,2-PB or cis-1,4-PB chains. As the SAFT model has been parametrized for experimental chain lengths (rather than the simulated chain lengths), this represents an excellent degree of transferability between molecular weights. In addition, from the slope of the lines in Figure 4, the coarse-grained models for squalane and PB correctly capture the isobaric thermal expansion behavior, a good indicator of transferability of nonbonded bonded potentials. We note that thermal expansion is often represented poorly within bottom-up structure-based coarse-grained models (particularly if these are obtained from iterative Boltzmann inversion 24 ) because of the strong density dependence of IBI nonbonded potentials.
We note in passing that, by design, the coarse-grained simulation model does not include angular terms in the force field. For longer chains these are expected to improve the structural representation of the chains (e.g., by providing a better radius of gyration). In the work of Morgado et al., 78 for short perfluoroalkylalkanes, bond length and bond angle potentials were successfully added to intermolecular parameters obtained by using the SAFT-γ Mie EoS. However, for the polymers used in this work, this approach cannot be used without a reparametrization of the nonbonded interactions (otherwise bulk densities are effected by addition of bond angle potentials). A slightly different approach was adopted in the work of Walker et al., 79 who have developed a fused-sphere SAFT-γ Mie force field for poly(vinyl alcohol) and poly-(ethylene). This approach uses the shape factor parameter within the fused-sphere SAFT-γ Mie theory, together with iterative Boltzmann inversion from fully atomistic simulations, to derive effective bond stretch and bend parameters, in addition to developing empirical correlations between the shape factor and the standard Mie parameters. Although this approach requires several steps, good agreement was obtained between the atomistic and coarse-grained radius of gyration and end-to-end distance distribution functions. A valuable extension to the SAFT-γ Mie theory would be provided by the addition of angle terms, which would allow a SAFT-based topdown coarse-graining methodology to better capture structural data in addition to thermodynamic data for longer chain polymers. However, this is not easy to achieve within the current formulation of the theory.
Mixture Compatibility. The measurements of cloud point temperatures for squalane in polybutadiene and cis-polybutadiene are plotted in Figure 5 along with fits obtained from SAFT-γ. If standard combining rules for Mie potentials could be used without k ij and γ ij parameters, then SAFT-γ Mie would provide an easy way to predict directly multicomponent phase diagrams for polymer−polymer and polymer−oligomer systems without the need for additional experimental data on mixing. Unfortunately, setting k ij = γ ij = 0, and using only combining rules for all unlike interaction Mie parameters (ϵ ij , σ ij , λ r ij , and λ a ij ), SAFT-γ Mie predicts complete immiscibility for squalane−polymer binary mixtures of polybutadiene and   Macromolecules pubs.acs.org/Macromolecules Article cis-polybutadiene. For k ij = γ ij = 0, we attempted also to fit the mixture and pure component data at the same time using combining rules for (ϵ ij , σ ij , λ r ij , and λ a ij ) values and also attempted to fit the same data with λ a ij fixed to a value of 6 (the usual value for the attractive exponent). Neither of these approaches was successful. Keeping the SAFT-γ Mie parameters from pure component systems and adjusting k ij to fit the cloud point curves provides a model that predicts partial miscibility, as seen experimentally, but with lower miscibility at high temperatures (LCST behavior), not improved miscibility at higher temperatures (UCST behavior). However, the γ ij parameter can change the system from an LCST to UCST mixture, as both Mie potential exponents change the cloud point curve's slope and curvature around the critical region. Optimized simultaneously, k ij and γ ij give enough flexibility to obtain partially miscible UCST models for both systems, in good agreement with the experimental cloud point temperatures. The optimized Mie potential parameters for the mixtures are listed in Table 2.
The final binary SAFT-γ Mie models for the mixtures sq/PB and sq/cisPB capture the cloud point curves and its UCST nature well. The exact UCST could not be experimentally measured due to the fact that samples at concentration of squalane above 80% (w/w) were not homogeneous enough to properly determine the temperature of the cloud point. However, it can be inferred that the experiment UCST, based on the curvature of the cloud point curve, is at a slightly lower squalane weight fraction than predicted by the SAFT model.
It is well-known that deuterium labeling can cause a small shift to phase boundaries. Because of the limited amount of deuterated squalane available, it was not possible to test this fully by measuring a complete phase diagram. However, we did confirm that deuteration of squalane has little effect on the UCST phase boundary with PB, with values of 53.8 and 49.2 ± 1°C at 15 and 20 wt % PB in squalane, respectively.
Surface Enrichment/Partial Wetting: Neutron Reflectivity. Initially, we checked for the presence of a wetting layer of squalane at the surface of PB for a 50% (w/w) mixture ( Figure 6). The wetting layer shows some temperature dependence, with surface segregation decreasing with increas-ing temperature. In addition, increasing the temperature to 45°C from 25°C, and then cooling back to the lower temperature, recovers the lower temperature results, suggesting that the results are equilibrium data and do not show hysteresis or suffer from significant problems with evaporation.
In Figure 7, we look at the surface enrichment of squalane for PB and PBcis films of different bulk concentrations. As the concentration of squalane in the mixture increases for both systems, we see a transition from surface segregation (up to a monolayer of squalane) to a full wetting layer. For PB− squalane this occurs at a bulk concentration of ∼40% (w/w). The presence of a wetting layer of higher scattering length density than the silicon substrate is evident as a shift in critical edge of R(Q) to Q > 0.01 Å −1 . The thickness of the wetting layer continues to increase with increasing bulk concentration. (We note in passing that specular NR is only sensitive to a onedimensional concentration profile and not to any variation in concentration lateral to the surface.) For PBcis the transition to a wetting layer occurs at a lower bulk concentration of squalane (∼30% (w/w)). This result can be expected from the two binary phase diagrams where for fixed temperature the two-phase region extends to lower squalane concentrations for PBcis.
Within the two-phase region it might be expected that the thickness of the surface wetting layer of squalane would depend on the thickness of the bulk film. However, this  Macromolecules pubs.acs.org/Macromolecules Article appears not to be the case. In Figure 8 we show results for PBcis from films of thickness 45 and 145 nm and a thick film (>300 nm). In these measurements the surface layer appears to be largely independent of film thickness and so appears to be in equilibrium with a phase-separated bulk system. Surface Enrichment/Partial Wetting: Molecular Dynamics. While the full experimental PB chain lengths are not accessible to the two-bead per-monomer chains used in this study, MD simulations can still be used in a semiquantitative way to study surface enrichment. Initial MD simulations of the squalane−polybutadiene mixture were performed at 450 K using parameters directly obtained from the SAFT-γ Mie models (above) and the standard mixing rules (eqs 3−6, k ij = 0). As the concentration of squalane is increased, the mixture undergoes phase separation, as seen experimentally. At all concentrations, the presence of a polymer surface leads directly to migration of the squalane to the surface, as shown by the simulation snapshot in Figure 9. At lower concentrations of squalane, for miscible mixtures, simple surface enrichment occurs. This behavior changes to a wetting layer at higher compositions with the wetting transition occurring between 40% and 60% (w/w) of squalane. As discussed above, the surface enrichment of the lower molecular weight species can be an entropy-driven effect, arising from the greater number of chain ends present in the lower molecular weight species. This effect has been seen previously in lattice models and in neutron reflectivity studies of mixtures and also single species polymers with different molecular weight components. 12,13 We additionally performed a series of simulation sweeps at 450 K as a function of squalane concentration and mixing parameter k ij . Entropy-driven migration of squalane to the surface was found for all k ij values studied. Here, bulk compatibility occurs at all concentrations for k ij ≤ −0.01 and bulk immiscibility for k ij > +0.005. For miscible systems, equilibration times are very fast on the order of a few tens of nanoseconds. However, equilibration becomes much slower when bulk phase separation occurs (see the Supporting Information for time sequence snapshots). Phase separation typically occurred over 80 ns when started from an initially homogeneous system, but equilibration of the surface layer depends on slow diffusion of phase-separated regions of squalane to the surface through the polymer film (microseconds time scales). Consequently, the time required for equilibration also increases with increasing system size. A typical sequence of curves is shown in Figure 10 for an initial bulk squalane concentration of 40% (w/w) with k ij parameters ranging from −0.1 (most compatible) to +0.005 (least compatible). As the unlike interaction parameter, k ij , is reduced, we see a move from surface segregation to a full wetting layer, with the thickness of the surface enrichment layer influenced by the value of k ij . Thicknesses of between 2 and 7 nm are found for the 40% (w/w) system (see Figure 10). These are consistent with experiment and in many cases larger than molecular dimensions. For comparison, a fully stretched squalane chain is 4.8 nm long. (Concentration sweeps are included in the Supporting Information for k ij = 0.0, −0.00165, and −0.01 along with results for larger systems of 102K particles, including results for a full wetting layer at k ij = +0.02.) A useful semiquantitative comparison to the experimentally obtained NR surface profiles is provided when the energetic   . Squalane and polymer surface density profiles obtained by CGMD at 450 K for systems with differing unlike interaction parameter k ij , ranging from attractive (k ij = −0.1) to repulsive (k ij = +0.005). The partial density of oligomers is shown by bold lines, and the partial density of polymers is shown by dashed lines. Surface enrichment varies from a small concentration increase at the surface to a pure squalane surface layer. An initial constant global squalane concentration of 40% (w/w) was used for each simulation. Results for higher values of k ij are shown in the Supporting Information.
Macromolecules pubs.acs.org/Macromolecules Article unlike interaction is set to k ij = 0, as given directly by the combining rules. Here, the slightly higher temperature used for simulations to achieve fast equilibration and the lower molecular weight of the PB chains offsets the need for the correction factors k ij and γ ij required for the full SAFT-γ Mie mixture model. We show comparison plots in Figure 11 for a range of concentrations experimentally and from simulation. As the bulk concentrations within the simulations are strongly influenced by migration, we note that comparison should be made for curves with similar bulk volume fractions. For CGMD and experiment we see incomplete surface coverage at low concentrations, and with increasing bulk concentrations, we see the maximum surface volume fraction of squalane approaches unity before the surface layer thickens significantly. In both cases, at the highest squalane concentrations, the surface layer thickness exceeds the dimensions of a single squalane molecule, ∼3 nm. However, we note that the predicted profiles are slightly steeper than the experimentally determined profiles, at both low and high squalane compositions, which are likely attributable to chain conformations being slightly too collapsed (smaller radius of gyration for chains) within the SAFT-derived MD model and also noting we are simulating polymer chains of shorter length than experiment. As above, we note that potential improvements to the MD force fields could be made through the simulation of longer chains and also the addition of angle and dihedral potentials, which will slightly increase the radius of gyration of the chains. These changes are also expected to result in rougher interfaces.

■ CONCLUSION
SAFT-γ Mie models for the sq/PB and sq/cisPB mixture are presented, which capture both liquid melt densities over a large pressure and temperature range and the blend compatibility across a wide range of temperatures and compositions.
In addition, the SAFT-γ Mie theory was used to directly produce a coarse-grained molecular dynamics (CGMD) model for the pure PB and squalane components and PB/sq mixtures in the presence of free surfaces. The CGMD provided excellent predictions for the temperature dependence of densities for squalane and very good predictions for PB, showing good transferability over different state points. The PB/sq mixture simulations were used to provide insights into the oligomer surface concentration profiles measured with neutron reflectrometry. For the latter we saw surface enrichment at lower squalane concentrations and a full wetting layer at high concentrations. Simulations with the CGMD force field were able to capture comparable surface enrichment from small bulk concentrations to the formation of a wetting layer.
This work provides one of the first CGMD simulations of a partially miscible oligomer/polymer mixture developed using the SAFT-γ Mie equation of state. Using this approach, it should be possible to predict, at least semiquantitatively, the surface migration behavior of a wide range of oligomers and polymers from their pVT data and phase behavior alone. Limitations of the method and possible improvements to SAFT-γ Mie and the CGMD potentials were discussed.
Phase boundaries obtained for PB/d-sq using different heating and cooling rates; simulation time snapshots showing development of a surface wetting layer for a 102K bead CGMD simulation 40% w/w oligomer; squalane and polymer surface density profiles obtained by CGMD for a larger 102K bead system for different values of k ij and different bulk concentrations (PDF)

■ ACKNOWLEDGMENTS
We are grateful to the European Commission (through the Marie Sklodowska-Curie actions 606869-MICSED) for funding this project and to Procter and Gamble (Germany), our partner, for helpful discussions. We also thank EPSRC for Macromolecules pubs.acs.org/Macromolecules Article funding through grant EP/P007864/1. We thank STFC for provision of the ISIS neutron reflectometry facilities and beamtime award RB1510402. We thank Dr. Buddhapriya Chakrabarti (Sheffield University), Dr. Gabriela Schafer, Dr. Torsten Lindner, and Dr. Jan Claussen (P&G, Germany) for many helpful and stimulating discussions. The authors thank PSE Ltd. for the use of its numerical solvers provided by the gSAFT package.