Molecular Simulations and Mechanistic Analysis of the Effect of CO2 Sorption on Thermodynamics, Structure, and Local Dynamics of Molten Atactic Polystyrene

A simulation strategy encompassing different scales was applied to the systematic study of the effects of CO2 uptake on the properties of atactic polystyrene (aPS) melts. The analysis accounted for the influence of temperature between 450 and 550 K, polymer molecular weights (Mw) between 2100 and 31000 g/mol, and CO2 pressures up to 20 MPa on the volumetric, swelling, structural, and dynamic properties of the polymer as well as on the CO2 solubility and diffusivity by performing molecular dynamics (MD) simulations of the system in a fully atomistic representation. A hierarchical scheme was used for the generation of the higher Mw polymer systems, which consisted of equilibration at a coarse-grained level of representation through efficient connectivity-altering Monte Carlo simulations, and reverse-mapping back to the atomistic representation, obtaining the configurations used for subsequent MD simulations. Sorption isotherms and associated swelling effects were determined by using an iterative procedure that incorporated a series of MD simulations in the NPT ensemble and the Widom test particle insertion method, while CO2 diffusion coefficients were extracted from long MD runs in the NVE ensemble. Solubility and diffusivity compared favorably with experimental results and with predictions of the Sanchez–Lacombe equation of state, which was reparametrized to capture the Mw dependence of polymer properties with greater accuracy. Structural features of the polymer matrix were correctly reproduced by the simulations, and the effects of gas concentration and Mw on structure and local dynamics were thoroughly investigated. In the presence of CO2, a significant acceleration of the segmental dynamics of the polymer occurred, more pronouncedly at low Mw. The speed-up effect caused by the swelling agent was not limited to the chain ends but affected the whole chain in a similar fashion.


Potential Energy Parameters for Atactic Polystyrene (aPS) and CO2
An all-atom representation was adopted, using the potential energy parameters from the work of Müller-Plathe S1 for aPS, in conjunction with harmonic constants for bond stretching from the work of Ndoro et al., S2 while the EPM2 model was chosen for CO2. S3 Instead of a rigid representation, the fluctuations of the bond length of CO2 were kept as low as possible by using a harmonic potential with a high value for the constant kb. It was verified that in this way the bond length of CO2 deviated by less than 0.01% from its equilibrium value of 1.1490 Å.  (3 )

Sanchez-Lacombe Equation of State (EoS) Parameters Regression
In the lattice fluid representation used to derive the Sanchez-Lacombe EoS S4-S6 , each substance is univocally characterized by three macroscopic parameters * , * , * . The ) is written in terms of reduced variables , , , defined as follows ( , , represent the temperature, pressure and density of the system): The EoS is formally identical for pure components and mixtures, provided that the corresponding definition of the reduced variables (i.e. that for a pure component or that for the multicomponent case) is used: is the total number of components in the system. Mixing rules to obtain the macroscopic parameters and reduced variables of a mixture are the following: ∆ * expresses the characteristic binary interactions between species i and j and contains an adjustable binary parameter , to account for deviations from the geometric mean mixing rule: represents the volume fraction of component in closed packed conditions, and can be calculated from the knowledge of the mass fractions : is the number of lattice cells occupied by a molecule in a mixture. It is related to the corresponding value for the pure component, , and the close-packed molar volume of a lattice cell for pure , * , or for a mixture, * : The expression for the chemical potential to be used in phase equilibria calculations is given below: The pressure-volume-temperature data reported by Zoller and Walsh S7 for polystyrene was used as target to obtain a functional dependence of the SL EoS parameters on Mw. The comparison between experimental data and EoS results is shown in Figure S1 and the parameters obtained for each case are reported in Table S4. A good fit could be obtained at all conditions, especially at lower pressures. The maximum deviation between experimental data and EoS results was 1.1% at 110000 g/mol, 1.0% at 34500 g/mol and 9000 g/mol, and 2.5% at 910 g/mol. The average deviation was 0.4% in all cases.

S10
The Mw dependence of the parameters is shown in Figure S. The results for * and * were interpolated using eq (S14), while for * eq (S15) was used.
( ) = 1 + (S14) The values for the best fit constants are: = 750 ± 10 and = 0.063 ± 0.02 for * , = 470 ± 35, = 0.36 ± 0.04 and = 360 ± 20 for * , = 1.099 ± 0.003 and = 0.023 ± 0.002 for * . * is expressed in K, * in MPa, * in g/cm 3 and Mw in g/mol. In the high Mw limit, these relations yield the parameter set from Doghieri and Sarti. S9 Using these relations, the parameters corresponding to the Mw of the simulated systems were obtained and they are reported in Table 2 in the main text. interpolations obtained with eq (S14) for (a) * and (c) * and eq (S15) for * (b). Furthermore, in order to calculate the properties of aPS-CO2 systems, the value for the binary interaction parameter must be determined. Since this parameter is usually temperature dependent, its value as well as its temperature dependence were determined by the best fit of the sorption isotherms reported by Sato et al. S10 in the temperature range 373 -453 K. This data set was selected because the sorption isotherm at 453 K showed the best agreement with the simulated data. For CO2 the SL parameter set reported by Doghieri and Sarti S11 was used. The binary interaction parameter is assumed to be independent of Mw and concentration. A linear extrapolation of the parameter values to the temperatures used in the simulations yielded the following values of the aPS-CO2 binary interaction parameter: = -0.041 at 450 K, = -0.060 at 500 K, = -0.080 at 550 K. This temperature dependence, however, led to a deviation from the expected Arrhenius behavior of the solubility as a function of temperature in the high-temperature range, as can be seen in Figure S20. The expected trend could be obtained by extrapolating the values of using a quadratic temperature dependence instead of a linear one: = 1.69 • 10 − 1.84 • 10 + 4.44 • 10 (S16) In this way, a value of = -0.041 was used at 450 K, = -0.052 at 500 K, = -0.055 at 550 K.
The EoS predictions of volumetric properties are almost insensitive to this choice, while in the case of CO2 solubility, a maximum 15% difference is obtained at 550 K (where the difference between the two values of is the greatest). S12 Density Figure S3. Specific volume of aPS as a function of inverse Mw. Orange circles are simulations at 450 K, red diamonds at 500 K, brown squares at 550 K. Results of Fox et al. S12 at 490 K included for comparison. Figure S4. Effect of cooling rate on the polymer density at different Mw.

Thermal expansion coefficients
The temperature dependence of the density was assessed by calculating the thermal expansion coefficients at atmospheric pressure for the pure polymer systems using eq (S17). The results are reported in Table S5 and compared to the values calculated for all the experimental and simulated data sets reported in Figure 1 in the main text.
Temperature fluctuations provide a negligible contribution in the calculated uncertainty compared to volume fluctuations.

Radius of Gyration
Values for the root mean squared radius of gyration / were obtained for all the systems using eq (S18) and are reported in Figure S6 as a function of CO2 content.
In the previous relation is the position of the center of mass of a chain, the position of atom of mass along the same chain and is the total mass of the chain. Error bars in the plots represent the standard deviation of the mean value.
The Mw dependence of / is shown in Figure S7 by comparing the simulated results at 450 K with values determined by neutron scattering for monodisperse aPS ranging from 21000 to 1100000 g/mol at 393 K. S20 S16 Figure S6. Radius of gyration of the polymer chains as a function of CO2 concentration. Plots

Radial Distribution Functions between Polystyrene atoms
Experimental measurements of the RDF of carbon atoms in polystyrene showed peaks at 1.4-1.5, 2.5, 5, 6, and 10 Å. S21,S22 The first peak is intramolecular and stems out from correlations between next neighbors, both on the chain and on the ring. These two contributions are resolved in the radial distribution functions calculated for the simulated systems ( Figure S8): the peak at 1.39 Å represents the distance between carbon atoms on the ring, while the lower peak at ~1.5 Å includes the contributions coming from neighboring backbone carbons (bond length 1.53 Å) and from the ring carbon bonded to the backbone (bond length 1.51 Å). The second peak is also

S18
attributed to intramolecular correlations, in this case of second neighbors in the backbone and the ring. These two contributions are fused in the simulation results; however, it is possible to attribute the shoulder on the right side of the peak to the correlations originating from the backbone, as can be seen in Figure 5 in the main text, where the contributions from ring and backbone carbons have been calculated separately. Additionally, another peak at 2.8 Å can be recognized in the radial distribution functions of the simulated systems, which represents the correlations associated with carbons located on opposite sides of a phenyl ring. A peak is located at 3.8 Å, and it can be associated with third neighbors along the backbone. Additionally, a small feature is present at 3.2 Å, and it is associated with backbone carbon correlations. Figure S8. Radial distribution function of pairs of carbon atoms (Mw 2100 g/mol at 500 K).

S21
Radial Distribution Functions between CO2 and Polystyrene atoms Figure S11 and Figure S12 show the effect of temperature and CO2 concentration on the RDF of CO2-polystyrene atom couples. With increasing temperature, the height of the peaks systematically decreases, especially in the features at lower distance, as can be observed in Figure S11. On the other hand, increasing CO2 concentration does not have a discernible systematic effect ( Figure S12). Generally, in the higher Mw systems, peaks tended to be slightly higher and broader, but this was not a systematic trend. It should be mentioned that the results at the lowest concentration suffer from worse statistics than the other cases, therefore the curves are less smooth, and it becomes difficult to establish trends in some cases.  The symbols represent simulated data at 500 K. Molecular weight of 2100 g/mol is depicted in orange, 5200 g/mol in red and 31000 g/mol in brown. Solid lines show extrapolation to shorter and longer times with a mKWW function.   Figure S15. Segmental relaxation times for pure aPS as a function of temperature. Simulations are represented with circles: brown represents 31000 g/mol, red is 5200 g/mol, orange is 2100 g/mol. Lines are NMR measurements, S23 arbitrarily shifted to compare the temperature dependence. Blue corresponds to a sample of 10900 g/mol, yellow 2100 g/mol, red 1600 g/mol.  Table S8. Repeating units included in each subsection in the analysis of segmental dynamics as a function of the position in the chain. Subsection I is the chain end, X is the center. I II III  IV  V  VI  VII  VIII  IX Table S9. Relaxation times associated to the decorrelation of the Car -Har bonds located at different positions in the chain at 500 K. Subsection I is the chain end, X is the center.  Figure S16. Effect of CO2 concentration on the reorientational decorrelation of Car -Har bonds in different chain subsections, for two polymer Mw at 500 K. Dynamics of End-to-end Vectors Figure S18. Effect of CO2 concentration on end-to-end vector reorientation at 500 K and at different molecular weights. Yellow represents results for the pure polymer, red at 5.70 × 10 -3 gCO2/gpol, brown at 2.82 × 10 -2 gCO2/gpol, blue at 5.05 × 10 -2 gCO2/gpol. In subplot (d) purple

Mw
represents the 2100 g/mol system, light blue represents the 5200 g/mol system, green represents the 31000 g/mol system.

Henry's Law Constant
Widom's test particle insertion test was performed on the pure polymer systems and the calculated values for the excess chemical potential were used to evaluate the Henry's law constant with eq (S19): In Figure S19, shown. They were calculated using the following expression of the model for the infinite dilution solubility coefficient for a two components system. Based on the definition given in eq (S19), this corresponds to the inverse of the Henry's law constants obtained in the simulations.

Enthalpy of Sorption
The values calculated with the Sanchez-Lacombe EoS using with a linear temperature dependence at high temperature deviate from the expected linear trend, whereas, if the values calculated with eq (S16) are used, the trend is linear in the whole temperature range. This finding would suggest that the extrapolation of the values of to higher temperatures using a linear relation might be inaccurate, and therefore it was abandoned.

S34
Self-diffusion coefficients Table S10. Self-diffusion coefficients of CO2 at 450 K, 500 K and 550 K and of aPS at 550 K.
In the high Mw case, the polymer did not reach a Fickian diffusion regime in the time of the simulation, therefore the diffusivity could not be extracted.