Analysis of Viscoelestic Properties of Elastomers Using Molecular Modeling

The goal of this research is to quantify and analyze elastomer chain conformations and their role in affecting rubber tire viscoelastic properties through the development of novel numerical methods. The hypothesis is that changes in the statistical mechanics of chain conformations at the microscopic scale provide a direct molecular link towards quantifying macroscopic properties of elastomers. The molecular level investigation began with simulating chains in the unperturbed state to study the impact of chain size and temperature on overall chain size statistics. To understand changes in the chain conformations under stress, a numerical model was developed that encompasses effects from multiple forms of deformation. The probabilistic approach implemented in this work allowed for molecular level understanding of chain behavior. Flory’s Rotational Isomeric State (RIS) approach was used to generate numerous uncorrelated, isolated, random conformations of amorphous cisand trans-1,4polybutadiene single chains under unperturbed conditions of different molecular weights and over a range of temperatures. Probability density distributions of squared end-to-end distances of these chains were quantified to study size properties. Characteristic ratios were in good agreement with prior experimental and theoretical findings, and increased with increasing molecular weight of cis and trans chains, with this effect being more pronounced for trans than for cis chains. Chain swelling was observed on heating indicated by an increasing characteristic ratio with temperature and positive temperature coefficients for both cis and trans chains. Chain size and shape properties were mutually dependent, with most changes in shape occurring due to changes along the principal direction. A larger relative increase in probability density distribution of unlikely larger chains and a smaller relative decrease in probability density distribution of more likely smaller chains resulted in increased average chain size and characteristic ratios with increasing temperature. This has been termed the “taut conformation effect” and had a significant impact on chain swelling with heating. This effect motivated further work into investigating its presence in other polymers such as polypropylene and polystyrene which is discussed in chapter 4. After the unperturbed state of chain ensembles were analyzed through size and shape studies, the perturbed or deformed state of the ensembles were explored through molecular modeling techniques to exert and quantify external stresses. Uniaxial, equibiaxial deformation and shear were applied to unperturbed chain ensembles which resulted in changes in their probability density distributions and elastic free energy. The approach for computing changes in elastic free energy involved developing a probability-based numerical method that can be applied across multiple forms of deformation. In order to determine the accuracy of the numerical method, it was initially applied to generated Gaussian chains and compared against known analytical equations. The numerical results and known analytical solutions of Gaussian chains were in excellent agreement and hence the numerical model was extended to computing elastic free energy change, force, and stress on RIS cisand trans-1,4-polybutadiene chains. Compression forces were much greater than tension forces. Equibiaxial and uniaxial stresses were equivalent in a single extension direction, and greater than shear. Forces and stresses increased with deformation and showed dependence on chain volume and temperature. Significant variation was observed in moduli with chain repeat unit size while only minor variations were observed with temperature. Fewer repeat unit size chains correspond to lower molecular weight between cross-links causing a more tightly cross-linked chain network. This resulted in greater moduli as compared to chains of more repeat units. A slight linear increase in moduli with temperature of chains of the same repeat unit size was observed. Young’s and shear moduli computed from the numerical model were in good agreement with experimental results. The novelty in this approach is the ability to incorporate polymer-polymer and polymer-filler interactions. The “taut conformation effect” was identified as a significant contributor to chain size behavior for polybutadiene, and its presence in vinyl polymers, such as polypropylene and polystyrene, were investigated. Random conformations of numerous single chains of amorphous, isotactic polypropylene and polystyrene were generated using Flory’s RIS approach. Characteristic ratios were in decent agreement with prior experimental and theoretical results. These ratios increased with molecular weight and were higher for polystyrene than polypropylene. Chain heating resulted in shrinkage for both polypropylene and polystyrene, indicated by decreasing characteristic ratios and negative temperature coefficients, which were in decent agreement with experimental results. Probability density distribution and chain size subset analysis indicated that only the least probable long size chains or taut conformations showed a decrease in probability density distribution and characteristic ratio with temperature. The most probable medium size chains hardly showed any change in their probability density distributions and characteristic ratios with temperature, while short size chains showed marginal increase. Hence, much like polybutadiene chains, less probable taut chain conformations had a significant impact on the average chain size indicating “taut conformation effect”. The primary focus of this dissertation was analyzing chain conformation statistics. In addition to that work, rheology experiments were performed on asphalt systems, and models were developed in order to predict the experimental results. Rheological studies were done to analyze viscoelastic properties of asphalt binders from different sources and under various aging conditions. Time-temperature superposition (TTS) was applied to frequency sweep data to produce master curves, and discrete and continuous Maxwell models were applied to predict stress relaxation modulus, relaxation and retardation time distributions, and creep compliance. Moduli increased upon aging indicating increased binder stiffness. Zero-shear viscosities decreased with increasing temperature for all binders due to increasing molecular motion and flexibility. At shorter times, stress relaxation moduli and creep compliance were similar for all binders, but with increasing time, unaged binders relaxed more rapidly than aged binders. Low creep compliance at shorter times corresponded to the absence of any configurational re-arrangements of asphalt binders. TTS allowed for computation of creep compliance at very low temperatures (-18◦C) to predict Pressure Aging Vessel (PAV) aged binder stiffness and resistance to thermal cracking at such low temperatures. Two binder samples were analyzed and showed similar low temperature flexibilities. The low temperature stiffness results predicted by our rheological model were in good agreement with bending beam rheometer experimental results, hence corroborating the efficiency of our model. This work encompasses three significant contributions to the scientific community. For polybutadiene, polypropylene and polystyrene, the impact of less likely, taut conformations on the overall chain size was identified and termed the “taut conformation effect”. A probability-based, reliable numerical tool was developed to predict mechanical properties of elastomers subjected to multiaxial deformations. Finally, an example of applying rheology models to experimental data was shown by accurately predicting asphalt binder stiffness and resistance to low temperature thermal cracking.

The molecular level investigation began with simulating chains in the unperturbed state to study the impact of chain size and temperature on overall chain size statistics. To understand changes in the chain conformations under stress, a numerical model was developed that encompasses effects from multiple forms of deformation.
The probabilistic approach implemented in this work allowed for molecular level understanding of chain behavior.
Flory's Rotational Isomeric State (RIS) approach was used to generate numerous uncorrelated, isolated, random conformations of amorphous cis-and trans-1,4polybutadiene single chains under unperturbed conditions of different molecular weights and over a range of temperatures. Probability density distributions of squared end-to-end distances of these chains were quantified to study size properties. Characteristic ratios were in good agreement with prior experimental and theoretical findings, and increased with increasing molecular weight of cis and trans chains, with this effect being more pronounced for trans than for cis chains.
Chain swelling was observed on heating indicated by an increasing characteristic ratio with temperature and positive temperature coefficients for both cis and trans chains. Chain size and shape properties were mutually dependent, with most changes in shape occurring due to changes along the principal direction. A larger relative increase in probability density distribution of unlikely larger chains and a smaller relative decrease in probability density distribution of more likely smaller chains resulted in increased average chain size and characteristic ratios with increasing temperature. This has been termed the "taut conformation effect" and had a significant impact on chain swelling with heating. This effect motivated further work into investigating its presence in other polymers such as polypropylene and polystyrene which is discussed in chapter 4.
After the unperturbed state of chain ensembles were analyzed through size and shape studies, the perturbed or deformed state of the ensembles were explored through molecular modeling techniques to exert and quantify external stresses.
Uniaxial, equibiaxial deformation and shear were applied to unperturbed chain ensembles which resulted in changes in their probability density distributions and elastic free energy. The approach for computing changes in elastic free energy involved developing a probability-based numerical method that can be applied across multiple forms of deformation. In order to determine the accuracy of the numerical method, it was initially applied to generated Gaussian chains and compared against known analytical equations. The numerical results and known analytical solutions of Gaussian chains were in excellent agreement and hence the numerical model was extended to computing elastic free energy change, force, and stress on RIS cis-and trans-1,4-polybutadiene chains. Compression forces were much greater than tension forces. Equibiaxial and uniaxial stresses were equivalent in a single extension direction, and greater than shear. Forces and stresses increased with deformation and showed dependence on chain volume and temperature. Significant variation was observed in moduli with chain repeat unit size while only minor variations were observed with temperature. Fewer repeat unit size chains correspond to lower molecular weight between cross-links causing a more tightly cross-linked chain network. This resulted in greater moduli as compared to chains of more repeat units. A slight linear increase in moduli with temperature of chains of the same repeat unit size was observed. Young's and shear moduli computed from the numerical model were in good agreement with experimental results. The novelty in this approach is the ability to incorporate polymer-polymer and polymer-filler interactions.
The "taut conformation effect" was identified as a significant contributor to chain size behavior for polybutadiene, and its presence in vinyl polymers, such as polypropylene and polystyrene, were investigated. Random conformations of numerous single chains of amorphous, isotactic polypropylene and polystyrene were generated using Flory's RIS approach. Characteristic ratios were in decent agreement with prior experimental and theoretical results. These ratios increased with molecular weight and were higher for polystyrene than polypropylene. Chain heating resulted in shrinkage for both polypropylene and polystyrene, indicated by decreasing characteristic ratios and negative temperature coefficients, which were in decent agreement with experimental results. Probability density distribution and chain size subset analysis indicated that only the least probable long size chains or taut conformations showed a decrease in probability density distribution and characteristic ratio with temperature. The most probable medium size chains hardly showed any change in their probability density distributions and characteristic ratios with temperature, while short size chains showed marginal increase. Hence, much like polybutadiene chains, less probable taut chain conformations had a sig-nificant impact on the average chain size indicating "taut conformation effect".
The primary focus of this dissertation was analyzing chain conformation statistics. In addition to that work, rheology experiments were performed on asphalt systems, and models were developed in order to predict the experimental results.
Rheological studies were done to analyze viscoelastic properties of asphalt binders from different sources and under various aging conditions. Time-temperature superposition (TTS) was applied to frequency sweep data to produce master curves, and discrete and continuous Maxwell models were applied to predict stress relaxation modulus, relaxation and retardation time distributions, and creep compliance. Moduli increased upon aging indicating increased binder stiffness. Zero-shear viscosities decreased with increasing temperature for all binders due to increasing molecular motion and flexibility. At shorter times, stress relaxation moduli and creep compliance were similar for all binders, but with increasing time, unaged binders relaxed more rapidly than aged binders. Low creep compliance at shorter times corresponded to the absence of any configurational re-arrangements of asphalt binders. TTS allowed for computation of creep compliance at very low temperatures (-18 • C) to predict Pressure Aging Vessel (PAV) aged binder stiffness and resistance to thermal cracking at such low temperatures. Two binder samples were analyzed and showed similar low temperature flexibilities. The low temperature stiffness results predicted by our rheological model were in good agreement with bending beam rheometer experimental results, hence corroborating the efficiency of our model. This work encompasses three significant contributions to the scientific community. For polybutadiene, polypropylene and polystyrene, the impact of less likely, taut conformations on the overall chain size was identified and termed the "taut conformation effect". A probability-based, reliable numerical tool was developed to predict mechanical properties of elastomers subjected to multiaxial deformations.
Finally, an example of applying rheology models to experimental data was shown by accurately predicting asphalt binder stiffness and resistance to low temperature thermal cracking.

ACKNOWLEDGMENTS
I would like to thank my advisor Dr. Michael Greenfield for his constant encouragement and valuable guidance throughout my Graduate School career. He has always encouraged me to strive for the best, to think outside the box, and push the boundaries of my knowledge. Our meetings ranged from an hour update to an evening marathon leaving both of us exhausted, but they were always fruitful and fulfilling. Thanks to him, I had the opportunity to travel all over the USA presenting our research, to network with other professional scientists, and work on projects in various aspects of my field. As I move beyond the realm of Graduate School, he has been instrumental in guiding me on a path towards a professional career. The experience I have gained from him is invaluable and will hold me in good stead for the future. I want to thank Brenda and Sheryl from the Chemical Engineering Department and Dean's Office respectively for their timely help whenever I needed them.
Eigenvalues of the radius of gyration matrix quantify chain shapes along principal directions. Averaged shape measures differ between cis and trans chains, and most shape changes arise along the longest principal direction. Joint correlations between chain size and shape show they are mutually dependent.

Introduction
Vehicle tires are an important application of rubber worldwide. Methods that improve fuel economy by decreasing rolling resistance without compromising wear resistance and traction have been reported [1,2,3,4]. Rolling resistance results from the energy a tire absorbs as it revolves and deflects when in contact with the road.
This paper takes steps toward understanding how elastomer chains contribute to rolling resistance. Its basis is an assumption that rolling resistance on the macroscale connects directly to energy losses from changes in chain conformations on the microscale. This paper is part of an overall project to quantify how differences in chain conformations in the presence of interactive fillers lead to entropy differences during deformation, which lead to changes in energy losses that con-tribute to rolling resistance. In this approach, studying and understanding the statistical mechanics of chain conformations provide a molecular link toward understanding the role of chain conformations in determining rolling resistance.
A rubber tire comprises one or more different types of elastomers such as styrene-butadiene rubber (SBR), polybutadiene, or polyisoprene. In addition to elastomers, tires include materials such as reinforcement fillers, curing agents, pro-  [9,10], and torsion angle probability distributions were similar within single chain and bulk structures. Kim, Misra, and Mattice [14] used that same methodology to study amorphous structures of trans-1,4-polybutadiene, isotactic 1,2-polybutadiene, and a random butadiene copolymer. They found few CH-CH 2 torsions in cis configurations. Gestoso et al. [15] used end-bridging Monte Carlo simulations to simulate cis-1,4-polybutadiene in the melt and as single chains.
They found similar end-to-end distances, single-chain structure factor, and torsion angle distributions for both cases.
Complementary results have been found for bulk systems of polybutadiene from molecular dynamics simulations. Predictions include density, dynamic structure factor, and torsion angle distribution [16,17,18,19,20,21,22,23]. Smith and Paul [24] developed a united atom (UA) force field for molecular dynamics simulation of 1,4-polybutadiene systems based on ab initio quantum chemistry.
They computed characteristic ratios for chains within bulk polybutadiene. Smith et al. [20] quantified differences among neutron spin echo scattering, molecular dynamics simulation, and the Rouse model, noting non-Gaussian distributions of segment-segment distances.
Some studies have focused on correlations in conformational changes. Gee and Boyd [18] found similar conformational transition rate activation energies in single chain and bulk systems, while torsion angle correlation rates showed a higher activation energy in the bulk. A cooperative kinetics approach models how relaxation rate about one torsion angle in polybutadiene impacts relaxation of neighboring torsions [25,26]. Correlations of this sort were found in a dynamic RIS model that was parameterized by MD results [27].
Experimental data may be interpreted to obtain chain conformation measures. Danusso et al. [28], Moraglio [29], and Abe and Fujita [30] measured intrinsic viscosity of cis-1,4-polybutadiene under theta conditions. Mark [9] extracted characteristic ratio results from these data. Mark [10] also cites intrinsic viscosity experiments published later [31] to determine the characteristic ratio of trans-1,4-polybutadiene. Hadjichristidis et al. [32] measured intrinsic viscosity under theta conditions to determine the characteristic ratio of a polybutadiene of 36% cis, 57% trans, 7% 1,2-vinyl composition. Gkourmpis and Mitchell [33] measured the structure factor of polybutadiene by neutron scattering and compared to predictions from the RIS model. Agreement between these methods was better at small wavevector (longer distances), which led them to alter the chain parameters in their model on the basis of improved agreement of predicted scattering results.
In this work we have studied size and shape properties of random chain conformations of polybutadiene by calculating the distribution of single chain conformations that arise under theta conditions. Numerous (> 10 5 ) uncorrelated random conformations of isolated cis-and trans-1,4-polybutadiene single chains were computed under unperturbed conditions. Using a single chain in each computation is justified because a flexible polymer surrounded by the same polymer takes on the same average shape as a single random polymer chain in a theta solvent or a melt [34], as confirmed in cis-1,4-polybutadiene simulations [13]. The polybutadiene chains were generated using the RIS model. Each chain realization in RIS provides an independent sample. While the standard molecular dynamics and Monte Carlo methods provide sequences of related states, the small changes that occur in each step lead to correlations that must be relaxed to sample an equilibrium distribution. The RIS method offers an advantage of generating a much larger number of uncorrelated random chain conformations in a computationally cheap manner. Using the ensemble of single-chain configurations, characteristic ratios were computed for cis-and trans-1,4 polybutadiene chains at different chain lengths and over a range of temperatures. Probability density distributions of the chains at different temperatures indicate a previously unknown result: relatively large increases in the small absolute probabilities of the most stretched chains are responsible for increases in average chain size (i.e. chain swelling) with increasing temperature. Chain shapes at different chain lengths and over a range of temperatures were also studied. Finally, joint probability correlations between chain size and shape for cis-and trans-1,4-polybutadiene chains were considered.

Methodology
In the RIS approximation, torsions about backbone bonds are treated as ex- At each condition, property averages used an equal weighting for each chain. This is appropriate because relative Boltzmann-weighted probabilities were taken into account while generating the chain conformations. Further calculations obtained correlations among chain size, length, and shape.

Chain Geometry and Internal Energy
The positions of all atoms in the chain are defined by the position of the chain start, the overall orientation of the chain, and the internal coordinates, i.e. the bond lengths, bond angles, and torsion angles. The chain start was always placed at (60, 60, 80)Å within a box of edge lengths (137, 137, 137)Å that used periodic boundary conditions. The angles that determine the overall chain orientation were chosen randomly. These choices for the chain start and orientation should not have any effect on the results. Bond lengths and bond angle supplements used in our computations were obtained from Mark [9,10] and are shown in Tables 1 and 2. Abe and Flory [11] used the same values in their calculations. For the C-C=C-C double bond, the torsional angle (φ) is zero (trans) or 180 • (cis).
The total energy of each system was calculated as a summation of the torsional energy [8] and the dispersion interactions between non-bonded atoms. Interactions among atoms separated by 3 or fewer bonds were included within the rotational isomeric state-dependent torsional energies. Atoms separated by three or more bonds contributed to the non-bonded interaction energy, which was computed using the Lennard-Jones (6-12) potential [35]. This short range interaction between widely spaced atoms was used to avoid direct chain overlap. Its range was restricted by using a quintic spline that started at r ij = 1.45σ ij and brought the interactions to LJ = 0 and d LJ /dr = 0 at r ij = 2.33σ ij . The scaling of average chain size with the molecular weight is checked below to confirm the chains remain unperturbed.
Every conformation of polybutadiene generated in our work has fixed bond lengths and bond angles, and thus the bond energies and angle energies are independent of conformation.

Chain Generation
Transformation matrices convert bond vectors from one reference system to another [8]. These orthogonal transformation matrices were used to determine atom positions within each single chain of polybutadiene from the internal coordinates. Four transformation matrices were used per repeat unit of polybutadiene: three for the C-C single bonds and one for the C=C double bond.
Each polybutadiene chain was built in an atom-by-atom manner by using the states chosen for three torsional angles per repeat unit (φ i , φ i+1 , φ i+3 ) around the three single C-C bonds, as defined in figure 1. Torsion angles φ i and φ i+3 about bonds i and i + 3 affect positions of the pendant hydrogen atoms (H i , H i , H i+3 , H i+3 ) attached to the backbone atoms C i and C i+3 . They also directly affect positions of the next atoms along the backbone (C i+1 , C i+4 ). Torsion angle φ i+1 directly affects positions of two carbon atoms (C i+2 , C i+3 ) and two hydrogen atoms (H i+1 , H i+2 ). The chain start and end were hydrogen atoms in place of C i−1 and C i+4 for the first and last repeat unit. Figure 1: Trans-1,4-polybutadiene structure showing bonds, bond angles and torsion angles. Atoms from C i to C i+3 and their pendant hydrogens comprise a single repeat unit. Numbering employs Flory's convention [8].
The rotational isomeric states (i.e. torsion angles) for each repeat unit were selected here by choosing randomly according to probabilities that account for relative energies of adding two consecutive torsion angles. Statistical weight matrices [8] that incorporate these relative energies for consecutive torsion angles were sug-gested by Mark [9,10] for 1,4-polybutadiene systems. The same set of matrices and statistical weights were used here, and chains were generated based on Boltzmannweighted probabilities. Baysal et al. [25] found that adjacent torsions show the most angle-angle correlation, with torsions further away correlated to some extent.
The partition function [8] incorporates each possible combination of rotational isomeric states of a chain. The pairwise probability of a single conformation equals its contribution to the partition function, divided by the partition function. Creating single chains in our simulations essentially creates realizations of the ensemble represented by this partition function. A linear congruential random number generator was used for the simulations. A different integer seed between 0 and 2 31 − 1 was used for each condition of chain length and temperature.

Chain size and shape parameters
An important chain size parameter is the squared end-to-end distance r 2 , which was calculated as where r x , r y , r z are the x, y and z coordinates of the end-to-end distance vector r.
The hydrogen atoms of the chain start and end defined the end-to-end distance, with periodic boundary conditions taken into account to ensure the proper image of each atom position was employed. The squared radius of gyration (r 2 g ) relates to the distance of each atom in the polymer chain from the center of mass, m j is the mass of atom j, r j = (x j , y j , z j ) T is the position vector of atom j of a polymer chain, r com = (x com , y com , z com ) T is the position vector of the center of mass of a polymer chain, and N is the total number of atoms in the chain. Proper periodic image calculations were included here as well.
Theodorou and Suter [36] used the eigenvalues (λ 1 , λ 2 , λ 3 ) of a radius of gyration matrix (S) to quantify contributions directed along the three principal directions (eigenvec-tors) of a chain conformation, where and other terms are defined analogously. The overbar indicates an average over all atoms in a single chain conformation. The radius of gyration matrix was transformed to a principal axis system, which diagonalised it such that its eigenvalues were in descending order (λ 1 ≥ λ 2 ≥ λ 3 ). Eigenvalue λ 1 corresponds to the principal direction with the longest dimension while λ 2 and λ 3 correspond to secondary directions. This approach effectively represents the size of a polymer chain using a rotational ellipsoid with a different size in each direction, rather than with a hollow sphere of radius r g that has the same mass and moment of inertia as the polymer chain. The squared radius of gyration equals the sum of the three eigenvalues, [36] The number-weighted squared radius of Eq. (5) differed by less than 1% from the mass-weighted squared radius of Eq. (2). Computing the radius of gyration matrix, Eq. (3), also enabled quantifying chain shape. The chain shape parameters studied were b/r 2 g (asphericity or deviation from spherical shape), c/r 2 g (acylindricity or deviation from cylindrical shape) and κ 2 (relative shape anisotropy), defined by Asphericity b/r 2 g equals zero when all dimensions are equal and goes to 1 when the principal direction is much larger, λ 1 λ 2 , λ 3 . Acylindricity c/r 2 g goes to zero when the secondary directions are equal lengths.

Chain size
The range of chain lengths allows for a limited test of the scaling of root-meansquared end-to-end distance. Figure 2 depicts the simulation results for r 2 1/2 at each chain length, using the number of repeat units rather than molecular weight or number of backbone bonds. Dashed lines indicate the best fits using a scal-  Characteristic ratio (C n ) [8,37] is defined as the ratio of mean squared endto-end distance of a real chain under the theta condition to that of a freely jointed chain with the same bond lengths and number of bonds. This condition is relevant to an elastomer because conformations in the bulk polymer above the glass transition resemble those in a theta solvent [34]. This assumption has been made by others within simulating polybutadiene chains to interpret neutron scattering data [33] and it was found to be accurate in a simulated melt [13]. Because two bond lengths arise in polybutadiene, they were incorporated separately as [10,11] C n = r 2 0 n s l 2 s + n d l 2   polybutadiene (open) using r 2 0 (•) and 6 r 2 g 0 ( ). Experimental results for cis-1,4-polybutadiene ( ) [9,28,29,30] , trans-1,4-polybutadiene ( ) [10], and mixed polybutadiene ( ) [32] are shown at 1/n = 0.001 . Prior model results for cis-1,4-( * ) [9,11,12,15,23,24], trans-1,4-(×) [10,11,12,24], and mixed polybutadiene ( ) [33] are shown at 1/n = 0.003. These symbols are used throughout unless otherwise specified. Figure 3 shows that calculated characteristic ratios are in good agreement with the experimental and prior computed characteristic ratios in the limit of infinite chain length for cis-1,4-polybutadiene chains. The only large simulation discrepancy is C ∞ = 3.76 from Kajiwara and Burchard [12]. They used a Monte Carlo simulation that sampled 10 3 to 10 4 of the rotational isomeric states for chains of up to 100 repeat units. Though they used the states proposed by Mark [9,10], their predicted characteristic ratio was about 25% lower than Mark's RIS result.
Monte Carlo calculations [15] using a full potential at 413 K led to C ∞ = 4.7, with a slightly lower value for continuous unperturbed chains; the latter are isolated chains that are subject to the full torsional potential, rather than divided into discrete rotational isomeric states. Molecular dynamics simulations [23] using the same potential and temperature led to C ∞ = 4.8. These show that the RIS parameters of Mark lead to slightly more extended chains than a full potential energy function, and temperature effects are expected to expand the difference to some extent. The lower value C ∞ = 4.6 for a mixed polybutadiene [33] was attributed by Gkourmpis and Mitchell to how their simulations incorporated torsional correlations. The lowest experimental value, C ∞ = 4.5, corresponds to a measurement by Abe and Fujita [30] at 332.85 K, which we interpret using a viscosity constant that is consistent with the average characteristic ratio of 4.9 that Mark reports [9] on the basis of other experimental data [28,29,30]. This value is also an outlier in figure 4.
For trans-1,4-polybutadiene chains, the calculated characteristic ratios are slightly higher than most prior computed values in the limit of infinite chain length.
The characteristic ratios increased with increasing chain length for both cis and trans chains. In both cases, chains of 120 repeat units were not yet at the long chain limit. This contrasts with an earlier demonstration that 30 to 40 repeat units were enough to reach the long chain limit in shorter Monte Carlo simulations [12].
The higher characteristic ratio for trans chains indicates a greater chain extension. This is potentially a consequence of the greater distance spanned between the carbon atoms bonded to the double bonded carbons, as noted by Mark [10] for a low energy state.
The characteristic ratio increased with temperature for both cis and trans chains, as shown in figure 4, and the increase was larger for trans than for cis polybutadiene chains. This indicates swelling of the average chain size upon heating. The trends are consistent with the prior modeling and experimental C ∞ results from the literature, particularly given that the n = 50 unit chain has a characteristic ratio somewhat below the long chain limit. The dependence on temperature may be compared to results from force measurements on swollen networks. Mark [9] reports d ln r 2 0 /dT = 0.40 × 10 −3 K −1 over 323 to 363 K for 1-4,cis-polybutadiene.
Studies that included temperatures of 298 to 338 K [31] confirmed an average of 0.41×10 −3 K −1 ; an earlier average [38] of 0.39×10 −3 K −1 was also cited. Crespi and Flisi [39] report slopes d ln r 2 0 /dT for cis chains that range from 0.12 × 10 −3 K −1 to 0.45 × 10 −3 K −1 over 293 to 333 K. Slopes increased with increasing cross link density, with decreasing chain extension, and with decreasing temperature. Mark [10] reports d ln r 2 0 /dT = −0.65 × 10 −3 K −1 for trans-1,4-polybutadiene but also notes experimental difficulties with the swollen network remaining stable, writing that "(these) results . . . must therefore be considered only qualitative estimates." Followup work [31] refined this uncertain average to −0.55 × 10 −3 K −1 . Gkourmpis and Mitchell [33] report a slight decrease in C ∞ with increasing temperature in their model calculations for a mixed polybutadiene. Our results may be compared by assuming n and l are independent of temperature, leading to Fits over 273-343 K for the 50-unit chain lead to 0.48×10 −3 and 1. In the limit of long chains without long branches, the mean squared radius of gyration r 2 g 0 should equal 1/6 of the mean squared end-to-end distance r 2 0 [8]. Figure 3 shows the ratio r 2 0 / r 2 g 0 was higher than 6 for shorter trans chains and decreased to 6 for longer chains. The ratio was slightly higher than 6 for cis chains at all chain lengths. Figure 4 shows that the ratio r 2 0 / r 2 g 0 was almost independent of temperature for cis chains, whereas for trans chains it increased with increase in temperature.
The probability density distribution of the squared end-to-end distance was calculated and compared with the probability density distribution for a Gaussian chain [8,37], The segments of each freely jointed chain in such an ensemble can be considered as performing a random walk in three dimensions with the only constraint being that each segment must be joined to its neighbors with a fixed bond length [8,37]. to 8000Å 2 , and even beyond these ranges the differences were small. We classify these as medium sizes within the entire range of the chain size distribution. These medium size chains have the highest probability of occurence. While smaller chains have a comparable probability density, they span a much smaller range of squared end-to-end distance. The Gaussian model predicted higher probability than the simulation results for shorter chains (size range of around 10 to 200Å 2 for cis and trans) as well as for longer chains (size range of around 3500 to 7000Å 2 for cis and around 8000 to 14000Å 2 for trans). Simulation results showed slightly higher probability than those predicted by the Gaussian model within parts of the medium size range. Cis chains showed a greater deviation from the Gaussian model compared to trans chains.
A similarity between conformations of RIS and Gaussian chains has been noted previously. Kajiwara and Burchard [12] calculated similar radii of gyration and hydrodynamic radii for RIS chains and Gaussian chains. Figure 5 shows that the similarities extend across the entire distribution of chain sizes. The temperature dependences of the probability density distributions of chain sizes for cis-and trans-1,4-polybutadiene are shown in figure 6. Squared end-toend distance has a much wider distribution than the squared radius of gyration.
Smaller size trans chains were slightly more probable at lower temperatures than at higher ones. Cis chains showed probabilities more independent of temperature. This temperature dependency of the chain size distribution appears to conflict with the temperature dependence of characteristic ratio. The characteristic ratios increased with temperature, as shown in figure 4, though figure 6 suggests similar probabilities for the most probable chains (i.e. low and medium size chains) at different temperatures. There was a larger relative increase, though smaller absolute change, in the small probability density for larger chain sizes with temperature, as compared to a smaller relative decrease in probability density for smaller chain sizes. While only small absolute changes in probability densities with temperature were observed for all chain sizes, the characteristic ratio showed a steady increase in average size with temperature. This effect was more pronounced for trans chains than for cis polybutadiene chains.
To examine this effect, characteristic ratios were calculated using only subsets of the chain size distribution, with results shown in figure 7. Chains with squared end-to-end distance r 2 < 600Å 2 were considered as smaller chains, those with squared end-to-end distance in the range of 600 to 3000Å 2 were considered as medium size chains, and those with r 2 > 3000Å 2 were considered as larger chains.
These squared distances correspond approximately to where the distributions from the RIS calculations and the Gaussian model cross each other in figure 5. A small decrease in characteristic ratio with increasing temperature was observed for smaller chains, while very slight rises in characteristic ratio were found for medium size chains. Characteristic ratios increased with increasing temperature for larger chains, and the increase was much more prominent for trans than for cis chains.
The fraction of chains that were within each range changed little with temperature, as shown by the inset in figure 7. Increases in characteristic ratio with temperature (polymer chain swelling, figure 4) can thus be attributed to the size increases of the relatively few extended and taut conformations, rather than expansion uniformly across conformations of all sizes. Despite there only being a small change in low probability conformations, these resulted in an increase in average chain size. We note that a sufficiently large sample size is required to obtain a sufficient number of these low-probability conformations that, by definition, occur rarely. The greater increase of characteristic ratio with temperature for larger chains, as shown in figures 4 and 7, indicates that this previously unreported "taut conformation effect" was more prominent for trans than for cis polybutadiene chains. This effect of rare extended chains is likely the cause of the difference in the temperature dependence of characteristic ratio in our results compared to those found in prior studies. Using a large enough sample set enabled larger numbers of rare extended conformations to be sampled.   repeat units. The inset shows the number fraction of (from top) medium, smaller, and larger chains.

Chain shape
Ensemble averages of chain shape parameters were obtained in order to quantify shape variations among polybutadiene chains. Since each chain establishes its own principal axes, the analysis uses a different coordinate system for each chain.
The results thus emphasize the deviations of each chain from a symmetric shape.
Rotation differences between the principal axes and the original (x,y,z) coordinates are not important and were not taken into account when combining the results into averages and distributions.
The eigenvalues λ 1 , λ 2 , and λ 3 of the radius of gyration matrix indicate the extents of orthogonal principal axes that span the region occupied by a chain in primary and secondary directions. Ratios of eigenvalues thus indicate if chains are being stretched or compressed. Figures 8 and 9 show the eigenvalue ratios as functions of inverse of chain length and temperature, respectively. These calculations were carried out at 343 K and for 50 repeat units respectively. Figure 8 shows that trans chains were more stretched than cis chains along the principal direction, though the distinction may vanish in the limit of long chains.
For trans chains, the extent of stretching decreased slightly with increasing chain length. For cis chains, the extent of stretching was larger with increasing chain length. The change in ratio between the two secondary directions followed the same trend as the principal direction but more subtly. These behaviors indicate that at the same chain length, trans chains were slightly less spherical than cis chains. Cis chains were slightly more elongated with increasing chain length, while trans chains were slightly less elongated. At long chain lengths, trans and cis chains have similar elongational deviations from a spherical shape. The changes in the ratio between the two secondary directions were small compared to the changes in average chain size r 2 g 0 . Figure 9 shows there was little or no variation in relative chain extent with temperature for cis chains. The principal direction ratio increased slightly with temperature for trans chains, while minor variations arose in ratio between the two secondary directions. This shows that as the temperature increased, trans chains were slightly more stretched along the principal direction and thus were slightly less spherical. Little or no variation in relative chain extent and thus in shape was found for cis chains as a function of temperature. Chain shape trends shown in  of chain length and temperature respectively. An asphericity factor (b/r 2 g ) of 0 suggests a spherical shape and 1 suggests a rod-like shape, while an acylindricity factor (c/r 2 g ) of 0 suggests a round cross section and 0.5 suggests a more flat cross section normal to the longest axis. κ 2 of 0 suggests a rod-like shape whereas 1 suggests structures of tetrahedral or higher symmetry [36]. length. This change in shape was more subtle for trans chains than cis. This behavior followed the same trend shown in figure 8 for the eigenvalue ratios. Figure 11 shows that cis chains exhibited little or no change in shape with temperature. Trans chains were slightly less spherical with increasing temperature.
This behavior followed the same trend shown in figure 9.
The relative shape anisotropy followed the same trend as asphericity as func-

Joint correlations in size and shape
Joint correlations between chain size and shape were studied to determine if their variations with chain length and temperature were independent or correlated properties. Cis and trans chains showed similar joint correlation behavior, with correlation and anti-correlation between chain size and shape occuring to a greater extent for trans chains as compared to cis chains. Visualizations that animate rotations of these three-dimensional plots are available as supplementary material. Differences P (b/r 2 g , r 2 g ) − P (b/r 2 g )P (r 2 g ) and P (c/r 2 g , r 2 g ) − P (c/r 2 g )P (r 2 g ) of 0 indicate that size and shape are completely in-dependent of each other, i.e. they act as mutually exclusive events. A positive difference indicates correlated events, while negative indicates anti-correlation.
For small rod-like chains, which arise less typically than average, figure 12 indicates some anti-correlation between size and shape. Small chains were nearer to spherical in shape, and high correlation between chain size and shape was observed for them. One physical interpretation is that for a chain of this length to have its ends near each other, its intermediate sections must expand in more htan one Cartesian coordinate direction. A hairpin turn that would enable conformations to emphasize a single direction is inconsistent with the allowed rotational isomeric states. For medium size chains, some correlation was found for chains that are near rod-like, while notable anti-correlation was found for more spherical chains. With the increased separation between chain ends at a fixed backbone length, the chain conformation apparently can emphasize a single direction more explicitly. Rodlike large chains showed correlation between chain size and shape, i.e. becoming more likely when a large distance between ends did not require a circuitous chain contour. Figure 13 shows size-shape correlations for acylindricity in cis chains. Small chains showed good correlation while being somewhat round in cross section.
Medium size chains showed correlation for chains that were more flattened in cross section. For medium size chains with round cross sections, the relationship between size and shape became anti-correlated. Large chains showed minor correlation between chain size and shape with being nearly round in cross section.
These correlations are consistent with the same explanations described regarding asphericity. In total, different size and shape probability density distributions were found for cis and trans chains over different chain lengths and across a range of temperatures. Probability densities are related to the work required to alter chain size and shape, and thus different probability densities for cis-and trans-1,4-polybutadiene conformations indicate different extents of work that must be done in order to alter chain size and shape. Quantifying this deformation work is the subject of ongoing research.

Conclusions
Ensemble averages and probability density distributions of sizes and shapes of cis-and trans-1,4-polybutadiene chains have been quantified for isolated single chains under undeformed theta conditions using a Rotational Isomeric State approach. Such conformations are considered to be representative for a chain in its own melt.
Characteristic ratios were larger with increasing chain length for both cis and trans chains, and these were in good agreement with experimental and prior computed values (cis-1,4-polybutadiene), and slightly higher than prior computed values (trans-1,4-polybutadiene). Characteristic ratios were higher for trans chains than for cis chains and indicate greater chain extension, which could be due to a greater distance spanned between the carbon atoms bonded to the double bonded carbons.
A Gaussian model predicted higher probability than simulation results at shorter and longer chain sizes for both cis and trans chains. Simulation results predicted a higher probability than the Gaussian model at certain regions of medium size chains for cis and trans chains while at other regions of medium size chains, simulation and Gaussian results were in agreement.
Characteristic ratios increased with increasing temperature for both cis and trans chains, with trans chains showing greater temperature dependence. The dependence for cis chains matched prior findings, while for trans chains the C ∞ (T ) were consistent with direct reports of characteristic ratio but opposite in sign to the negative d ln r 2 0 /dT that have been reported from force measurements on networks. Despite the rise in characteristic ratio, only small absolute changes in chain size probability densities with increasing temperature were calculated.
Smaller chain conformations showed a smaller relative decrease in probability density with increasing temperature, as compared to a larger relative increase for improbable larger chain conformations. In a subdivided chain distribution, smaller chains showed a small decrease in characteristic ratio with increasing temperature, medium size chains showed little or no variation in characteristic ratio, and larger chains showed more significant increases in characteristic ratio with increasing temperature. This accounted for an increase in average characteristic ratio of cis-and trans-polybutadiene chains with increasing temperature. This newly reported "taut conformation effect" was more pronounced for trans-than for cis-1,4-polybutadiene chains. It indicates that the increase in chain size originates from the fewer chains that are largest in size, rather than from increases in size among chains across the size distribution.
With increasing chain length, trans chains were less elongated while cis chains were more stretched along the principal direction. At the same chain length, trans chains were slightly less spherical than cis chains. Cis chains were slightly less spherical with increasing chain length, while trans chains were slightly more spherical. At long chain lengths, trans and cis chains reach similar shapes. The extent of stretching and compression was greater along the principal direction than the secondary directions. With increasing temperature, trans chains were slightly stretched along the principal direction whereas cis chains showed little or no change in shape. Thus trans chains were slightly less spherical with increasing temperature, while little or no variation in shape was computed for cis chains.
Variations of the largest eigenvalue λ 1 of the radius of gyration matrix have the most significant effects on chain shapes: most changes in shapes arose from changes along the longest principal direction.
At longer chain lengths, both cis and trans chains showed similar asphericity.
Little or no variation was computed in acylindricity for either cis or trans polybutadiene chains. Relative shape anisotropy followed the same trend as asphericity as functions of both chain length and temperature for cis and trans polybutadiene chains.
Joint correlation studies between chain size and shape showed that they are mutually dependent properties. For asphericity, rod-like small size and spherical medium size cis chains showed anti-correlation between chain size and shape.
Spherical small size, near rod-like medium and large size chains showed correlation between chain size and shape. For acylindricity, medium size chains of flattened cross section, and small and large size chains of round cross section showed correlation between chain size and shape. Round cross section medium size chains showed anti-correlation between chain size and shape. Trans chains showed similar behavior as cis chains with correlation and anti-correlation between chain size and shape occuring to a greater extent.
Cis-and trans-1,4-polybutadiene show different size and shape probability density distributions, which imply different amounts of deformation work to alter chain shape and size. Quantifying this deformation work and its implications for mechanical properties, viscoelastic properties, and rolling resistance are the subject of ongoing work.

Supplementary Material
Visualizations that animate rotations of the three-dimensional plots of the size-shape correlations calculated for cis-1,4-polybutadiene, figures 12 and 13, are available as supplementary material which can be found online.

Acknowledgments
We thank the Ford Motor Company University Research Program for funding this research. Keywords: Uniaxial deformation, numerical method, Gaussian analytical, polybutadiene chains, stress-strain, moduli

Introduction
Polybutadiene is an important commercial polymer with pertinent applications in the automobile industry. It is such an essential ingredient in rubber tires that nearly 70% of polybutadiene manufactured worldwide is used in the production of vehicle tires. Additional uses include manufacture of golf balls, toughened plastics, shoe soles, gaskets, shock absorbers and others [1,2].
As a tire rolls, the stresses on the tread exert strains through both elastic and viscous mechanisms. The latter results in energy dissipation leading to rolling resistance, which is the viscoelastic energy lost at low frequencies during flattening and re-rounding of tires. When a vehicle is in motion, the tire tread flattens against the road which results in the elastomer chains undergoing a change in their conformations [3,4]. This leads to affine deformation of the elastomer system which changes the number of ways the elastomer chains can be arranged. The change in arrangement of elastomer chains affects the elastic free energy of the system, which is logarithmically related to the probability distribution of the chain end-to-end vectors [5]. Chains and particles in the undeformed portion of the tire tread relax.
The original distribution of the chains is restored by random fluctuations after deformation, and this change in elastic free energy requires work which is dissipated as heat leading to rolling resistance. Thus studying the chain conformations and their changes under deformation is of utmost importance to understanding the viscous mechanisms of this rubber system.
Prior experimental and theoretical work have been done to study changes in the molecular structure of rubber systems. NMR generally measures the average chain deformation and re-orientation on a local scale [6,7,8,9] and SANS computes the average radius of gyration of chain conformations in a system [10]. Dubrović et al. [17] and Valić [18] studied molecular structure, dynamics, and segmental motion of rubber system molecules at lower deformation i.e. ratio λ = L/L 0 = 1.5. Dubrović et al. [17] studied the effect of uniaxial deformation on the microstructure of natural rubber during irradiation. They found that increasing the irradiation dose decreases the degree of swelling of natural rubber.
A deformation of 1.5 introduced structural changes in the rubber matrix. Valić [18] studied orientational motion of natural rubber chain segments cross-linked by γ-irradiation under fixed uniaxial deformation of 1.5. Direction of applied force determined the dynamic behavior of chain segments and resulted in permanent spatial orientation of rubber chain segments in the direction of the force.
Starkova and Aniskevich [19] carried out uniaxial tension tests on silica-filled SBR rubber in order to determine the incompressibility limit of the rubber and the corresponding Poisson's ratio. They found that the rubber system was incompressible up to a deformation of 4.5 and corresponded to a Poisson's ratio of 0.5.
Beyond a deformation of 4.5 the Poisson's ratio decreased to 0.4.
Mark [20,21] used the Rotational Isomeric State (RIS) approach proposed by Flory [22,23] to study configurational statistics of single cis-and trans-1,4polybutadiene chains under theta conditions. In the RIS approximation, torsions about single bonds are treated as existing in one or more discrete rotational states with each of these states chosen to coincide with a region of low potential energy.
The RIS approach has been used as an effective tool to study single chain properties of polybutadiene by several researchers including Abe and Flory [24], Kajiwara and Burchard [25,26], and also by us [27]. In our previous work [27], we studied size and shape properties of amorphous cis-and trans-1,4-polybutadiene systems.
Single chains of polybutadiene under unperturbed conditions were generated, and single chain size analysis and shape analysis were performed. Characteristic ratios of both cis-and trans-1,4-polybutadiene chains increased on heating indicating chain swelling. Probability density distributions of the chains over a range of temperatures showed a previously unreported effect: increase in chain size on heating originated from the least likely, taut conformations rather than uniformly across the size distribution. This was termed as the "taut conformation effect".
In this work, we propose a novel numerical method for simulating mechanical properties of elastomer chains subjected to uniaxial deformation. The numerical method arises from a probabilistic approach for evaluating the elastic free energy change of chain end-to-end vectors under deformation. From previous work [27], we have available distributions of end-to-end distances of cis-and trans-1,4polybutadiene chains under unperturbed conditions that were generated using the RIS approach. Here we supplement those with numerical sets of end-to-end distances that follow a Gaussian distribution. These are used to confirm that the numerical deformation methodology implemented here obtains the known analytical stress-strain relationship for Gaussian chains. Computations were then performed to predict the stress-strain relationship of unoriented RIS chains subjected to uniaxial deformation by using changes in chain conformation statistics. This paper is part of an overall project aimed towards studying elastomer chain conformations and their role in computing rubber tire viscoelastic properties.

Methodology
The intent of the new methodology is to obtain mechanical properties from conformational statistics in a numerical approach that can incorporate directly the effects of specific polymer-polymer and polymer-filler interactions, such as functionalized end-groups with silica filler interactions. The probability-based approach developed here enables sampling over many configurations without the cost of molecular dynamics simulations.
The approach arises from the elastic free energy change of single chains of elastomers under deformation in order to quantify their mechanical properties.
Elastic free energy A el of a single chain under deformation with an end-to-end vector r is related to the probability density distribution of the end-to-end vector of the chain [5,23,28] where P ( r) is the probability density distribution of end-to-end vectors of a single chain, c(T ) is only a function of temperature T and k b is the Boltzmann constant. The change in elastic free energy of a single chain from an undeformed to a deformed state is The change in elastic free energy of an ensemble containing N chains was obtained by integrating over all chains. This was achieved by multiplying equation 13 by the probability density distribution of r over all r and integrating over all possible end group positions, In order to accurately compute the change in elastic free energy for chains having cross-links, it is necessary to distribute a front factor [29] or the ratio of mean- where (∆r) 3 = (∆r 0 ) 3 is the volume of a voxel, r 2 0 and r 2 represent undeformed and deformed chain end-to-end distances respectively, and r 2 0 is the average of the mean squared end-to-end distance of the undeformed chains. The P ( r) of chains were computed based on their voxel location along the x, y, z coordinate directions P ( r) = number of chains in voxel xyz at r volume of voxel xyz * N Deformation applied to elastomers induces strain on the systems, which leads to changes in elastic free energy. This induced strain leads to quantifiable forces and stresses. From thermodynamics, force is directly related to the change in elastic free energy with system dimension [29] and is given as where deformation or stretch ratio λ k = r k /r k,0 [30] applied in the kth direction leads to force f k in the kth direction. r k and r k,0 are deformed and undeformed chain end-to-end vectors respectively and L is V 1/3 where V is the volume of a single chain. The tensile stress acting on the system due to deformation is obtained by dividing the force with the cross-sectional area (A m ) whose unit vector is aligned with the applied force The choice of cross-sectional area, as actual or initial, results in either true or engineering stress [29].
Uniaxial deformation was applied to chain ensembles such that extension was along the x direction and compression was along the y and z directions. Affine deformation was assumed, i.e. deformation applied macroscopically is transferred microscopically and uniformly to every chain making up the ensemble [29,31].
Since polybutadiene is an elastomer, it was assumed to be incompressible with a Poisson's ratio (ν) of 0.5 [19,30]. No change in volume of the system was expected as a result of deformation thus λ x λ y λ z = 1. The deformations applied along the y and z directions are compressions such that λ y = λ z = 1/ √ λ x .
In this work true stress has been computed and cross-sectional area is represented by L 2 /λ k . Thus equation 18 can be written as The volume of a single chain (V ) is a function of its molecular weight and temperature, and was calculated using polybutadiene densities. Densities of cis-and trans-1,4-polybutadiene at different temperatures were obtained from prior experimental and theoretical results [32,33,34].
Tension force in the x direction can be obtained numerically from equations 15 and 17 as where λ is changing in all three directions i.e. extension along the x direction and compression along the y and z directions. For compression, where f y,z are compression forces in the y and z directions respectively.
Numerical equations for tensile and compressive stresses are obtained by dividing tension and compression forces with their respective cross-sectional areas, which are L y L z = L 2 /λ for tension and L x L y = L x L z = L 2 √ λ for compression Unperturbed RIS chain end-to-end distances (r 0 ) from previous work [27] were used as inputs for computing end-to-end vectors of chains before and after deformation. In order to evaluate equations 20 to 23, the end-to-end vectors had to be distributed over a three dimensional space. The chain conformations were randomly oriented over the (+x, +y, +z) octant of a sphere prior to imposing deformation. Distributing the chains over a single octant of a sphere assumes that the same distribution occurs along the remaining sections of the sphere, and it provides better statistics while requiring less computation time.
Deformation on the chains is implemented by applying the deformation ratio, representative of either extension (λ) or compression (1/ √ λ), to each of the vector components in their respective x, y, and z directions. First we considered randomly distributed unperturbed end-to-end vectors of the chains in spherical coordinates, (r x,0 , r y,0 , r z,0 ) = r 0 (sin θ cos φ , sin θ sin φ , cos θ) Each elevation or latitudinal angle θ [0,π/2] was chosen uniformly from a sinu-soidal distribution between [0,1] and longitudinal angle φ was chosen uniformly within [0,π/2). Probability density distributions for the unperturbed end-to-end vectors were obtained by averaging over all r 0 . After deformation, the end-to-end vectors of the chains were represented by (r x , r y , r z ) = r 0 (λ x sin θ cos φ , Probability density distributions for the deformed end-to-end vectors were obtained by averaging over the number of chains ending in the region around each r.
In order to verify the accuracy of the proposed numerical method, it was applied to Gaussian chains and the results were compared with analytical solutions.
The probability density distribution for the end-to-end vector r 0 of a Gaussian chain is given as [23] P ( r 0 ) = 3 2πC n (n s l 2 where r 0 is the end-to-end distance of a chain, C n is the characteristic ratio [22], n i is the number of backbone bonds of each type along a polymer chain, l i is the bond length, and subscripts s and d indicate single and double bonds where l s = 1.53Å and l d = 1.10Å for polybutadiene [20]. Gaussian chains were generated by repeatedly determining the distance r 0 that satisfies where ξ was chosen uniformly within [0,1). A linear congruential random number generator was used to generate each ξ with a different seed used in each system to ensure independent statistics.
The change in elastic free energy for an ensemble of N Gaussian chains is obtained by substituting equation 26 in equation 13 and averaging over N chains [5] ∆A el The mean squared end-to-end distances can be represented in terms of end-to-end vectors as Assuming the elastomer system to be isotropic in a state of rest [5] r 2 x 0 = r 2 y 0 = r 2 z 0 = The average mean squared end-to-end distances can be written in terms of deformation ratio as where i = x, y, z. The free energy change for Gaussian chains can be written in terms of deformation ratios [5] by substituting equations 29, 30, 31, and 32 in Under uniaxial deformation in the x direction, equation 33 can be written as [30] ∆A el and compared with results obtained from using equation 15 for Gaussian chains to determine accuracy of the numerical method for elastic free energy change computations.
Analytical tension (f x ) and compression forces (f y,z ) on the Gaussian chain ensembles were obtained from equations 17 and 33 and compared with equations 20 and 21 respectively to verify the numerical method accuracy for force computations.
Analytic tensile (σ xx ) and compressive stresses (σ yy,zz ) were obtained for the Gaussian chain model by dividing the tension and compression forces by their respective cross-sectional areas, leading to

Results and Discussion
A significant step towards determining the accuracy of the numerical method was comparing our numerical results with Gaussian analytical equations. In order to bring about such comparisons, Gaussian chains were generated using equation 27. Figure 14 shows good agreement in the probability density distribution of generated Gaussian chains with the Gaussian analytical probability density distribution.       This indicated that stiffness of the chains increased up on heating. These moduli were in good agreement with experimental results [38,39]. Textbooks on rubber elasticity [5,40] and polymer engineering [30] have reported Young's modulus of polybutadiene to be ∼1 to 3 MPa and our numerical results fall within that range.

Conclusions
The

Acknowledgments
We thank the Ford Motor Company University Research Program for funding this research.

Methodology
Part 1 focused on uniaxial deformation of polybutadiene chains. There, we developed a probability-based numerical method for computing the change in elastic free energy and stress-strain behavior of chains. In this work, we extend that methodology to study mechanical properties of chains subject to equibiaxial deformation and shear. Prior to being applied to non-Gaussian RIS chains, the accuracy of the method for these new deformation directions was tested against known analytical equations by using simulated Gaussian chains.
The probability-based numerical method being described here for the change in elastic free energy of an ensemble of N chains due to deformation is given as where (∆r) 3 = (∆r 0 ) 3 is the volume of a voxel, r 2 0 and r 2 represent undeformed and deformed chain end-to-end distances respectively, and r 2 0 is the average of the mean squared end-to-end distance of the undeformed chains. Each voxel represents a possible position of an end-to-end vector in three dimensional space. The P ( r) of chains were computed based on their voxel location along the x, y, z coordinate directions P ( r) = number of chains in voxel xyz at r volume of voxel xyz * N From thermodynamics, force is directly related to the change in elastic free energy with system dimension [12] as a result of deformation and is given as where deformation or stretch ratio λ k = r k /r k,0 [1] applied in the kth direction leads to force f k in the kth direction. r k and r k,0 are deformed and undeformed chain end-to-end vectors respectively and L is V 1/3 where V is the volume of a single chain. Stresses acting on the system due to deformation are obtained by dividing the forces with the appropriate cross-sectional areas (A m ) Biaxial deformation involves deforming the chains in tension along the x and y directions with compression along the z direction. For shear, we assume chains to be sheared along the x and z directions and unchanged along the y direction. Affine deformation was assumed [12,13] and since polybutadiene is an elastomer, it was assumed to be incompressible with a Poisson's ratio (ν) of 0.5 [1,14]. No change in volume was expected as a result of deformation thus λ x λ y λ z = 1. For biaxial deformation, the tension along x and y directions were represented by deformation ratios λ x and λ y respectively. Compression occurs along the z direction such that λ z = 1/λ x λ y . Deformation ratio for shear along the x direction is λ x or λ and along the z direction is 1/λ, with no change in the y direction i.e. λ y = 1.
For tensile stresses, the cross-sectional areas of interest are those whose unit vector is aligned with the applied force, and for shear stresses the cross-sectional areas are normal to the applied force. The cross-sectional areas for biaxial tension and shear For biaxial tension, substituting area from equation 46 or 47 into 45 leads to Similarly, substituting equation 49 into equation 45 for shear leads to The volume of a single chain (V ) is a function of its molecular weight and temperature, and was calculated using polybutadiene densities. Densities of cis-and trans-1,4-polybutadiene at different temperatures were obtained from prior experimental and theoretical results [15,16,17].
Forces can be obtained numerically by applying equation 44 to equation 42 for biaxial tension and shear Shear force in the z direction can be given as Stresses were numerically obtained by substituting equation 42 in equations 50 and 51 respectively for biaxial tension and shear The numerical method was tested by applying it to a distribution of Gaussian chains, for which analytical force and stress equations are available. Gaussian chains available from previous work on uniaxial deformation were oriented randomly over the (+x, +y, +z) octant of a sphere using combinations of θ and φ (part 1). The end-to-end vectors of the chains (r x,0 , r y,0 , r z,0 ) in spherical coordinates under undeformed conditions were (r x,0 , r y,0 , r z,0 ) = r 0 (sin θ cos φ , sin θ sin φ , cos θ) and after deformation (r x , r y , r z ) biaxial = r 0 (λ x sin θ cos φ , λ y sin θ sin φ , 1 λ x λ y cos θ) (r x , r y , r z ) shear = r 0 (λ sin θ cos φ , sin θ sin φ , The change in elastic free energy of an ensemble of N Gaussian chains is written analytically as [2] ∆A el Under biaxial deformation and shear, equation 64 can be written as and compared with results obtained from using equation 42 for Gaussian chains to determine accuracy of the numerical method for elastic free energy change computations.
Analytical forces on the Gaussian chain ensembles were obtained from equa-tions 44 and 65 for biaxial tension and compression Similarly, analytical shear force on Gaussian ensembles in the x and z directions can be obtained from equations 44 and 66 as The analytical force equations were compared with their numerical counterparts (equations 52, 53, 54, 55, and 56) to verify the accuracy of the force calculations within the numerical method.
Analytic stresses on Gaussian chain ensembles were obtained by dividing forces with the respective cross-sectional areas. This led to equations for biaxial tension and shear repeat unit size of 50 [11]. Following the numerical and analytical comparisons of Gaussian chains, the numerical method was extended to studying stress-strain behavior of non-Gaussian RIS chains.
Slopes of the linear regime of the respective stress-strain plots were used to compute Young's modulus (E) and shear modulus (G) [1,12]. The magnitude of the deformations imposed here (λ = 1 to 3.5) emphasized the need for utilizing finite strain ( ,γ) which can be related to engineering strain ( e ,γ e ) as [ For shear, engineering strain (γ e ) is related to λ as [1] γ e = λ − 1 λ (79) thus finite shear strain (γ) in terms of λ is Shear and Young's modulus can be related as [1,12] For elastomer systems of ν = 0.5, equation 81 can be simplified to G = E/3.

Results and Discussion
In order to test the accuracy of the numerical method, comparisons were done with analytical solutions for Gaussian chains under equibiaxial tension and shear. Figure

Conclusions
The aim of this work was to apply our new probability-based numerical method toward computing changes in elastic free energy and mechanical properties of elastomer systems under equibiaxial deformation and shear. Excellent agreement was observed between the numerical method applied to Gaussian chains and analytical functions, hence proving its correct implementation. This method was then applied to RIS generated cis-and trans-1,4-polybutadiene chains to predict their stress-strain behavior. Studies were done on chains of varying repeat unit size and temperature, and were compared with previous findings under uniaxial deformation (part 1). Elastic free energy change, force and stress increased with deformation for all three types of deformation. The tension force under equibiaxial and uniaxial deformation, and shear in the x direction were equivalent.
Compression forces were much larger than extension or shear forces. Force, stress, and ultimately moduli showed dependence on volume or temperature for chains of the same temperature or repeat unit size respectively. Similar tensile stresses were observed for chains under uniaxial and equibiaxial extension per direction of stretching, and were greater than shear stresses. Fewer repeat unit chains correspond to a much more tightly cross-linked network than more repeat unit chains, resulting in greater moduli. Good agreement in moduli was found between our numerical method and experimental results, with shear moduli being ∼1/3 of Young's moduli. Our numerical method in combination with the RIS method [22] has the potential to incorporate specific polymer-polymer and polymer-filler interactions, that cannot be accommodated easily for analytical models of chain configuration.

Acknowledgments
We  propylene and styrene respectively [2]. Their repeat unit structure can be writ-ten as -(CH 2 CHR) x -where x is the number of repeat units and -R is -CH 3 for polypropylene and -C 6 H 5 for polystyrene. The stereochemical relationship of consecutive -R groups along the polymer backbone determines the tacticity of the polymer. Isotactic and syndiotactic chains are stereoregular which tend to crystallize easily in bulk systems while stereoirregular atactic chains do not crystallize as easily, hence we have focused our studies on atactic chains in amorphous state.
Polymer chain conformations enable the computation of chain size properties such as characteristic ratio (C n ) [3,4] and its temperature coefficient [3]. Characteristic ratio can be defined as the ratio of mean-squared end-to-end distance of a real chain under theta conditions to that of a freely jointed chain with the same number of bonds and bond length. Temperature coefficient shows variation in chain size with temperature and can be computed from characteristic ratio.
The most common and cost effective method for experimentally determining average polymer chain sizes is through intrinsic viscosity experiments. Various molecular weight fractions of a polymer are dissolved in theta solvents and analyzed to determine intrinsic viscocity, which ultimately lead to mean-squared end-to-end distances of unperturbed chains [3,5]. Kinsinger and Hughes [6,7], Inagaki et al. [8], Nakjima and Saijyo [9], and Heatley et al. [10] performed such experiments on polypropylene under theta conditions. Using their data, Suter and Flory [11] computed characteristic ratios and temperature coefficients. Similarly, Krigbaum and Flory [12], Orofino and Mickey [13], and Altares et al. [14] carried out such experiments for polystyrene and reported intrinsic viscosities and mean-squared end-to-end distances. Flory [3] computed characteristic ratio from those data.
Orofino et al. [15] analyzed a sample of cross-linked bulk polystyrene over 393. 15-448.15 K. At each temperature, the sample was stretched to a desired elongation and allowed to relax until the applied tension appeared constant. Meansquared end-to-end distance and temperature coefficient were calculated from tension. Additionally, they performed intrinsic viscosity experiments on polystyrene in different theta solvents, at different temperatures and reported a similar temperature coefficient. Cotton et al. [16] performed neutron scattering experiments on amorphous polystyrene in a bulk system and determined mean-squared radius of gyration. Measurements were made for eight monodisperse polystyrene samples over a wide range of molecular weights. They compared the conformations in bulk to those in a theta solvent, and inferred that the chain dimensions in bulk were the same as those in a theta solvent. They also reported negative temperature coefficient for polystyrene.
Theoretical methods enable generating polymer chain conformations and studying their properties. The Rotational Isomeric State approach (RIS) proposed by Flory [3] has been a powerful theoretical tool for such investigations.
The RIS method is efficient for sampling single-chain conformations, as each chain realization in RIS provides an independent sample. While the standard molecular dynamics and Monte Carlo methods provide sequences of related states, the small changes that occur in each step lead to correlations that must be relaxed to sample an equilibrium distribution. The RIS method offers an advantage of generating a much larger number of uncorrelated random chain conformations in a computationally cheap manner [17]. In the RIS approximation, torsions about backbone bonds are treated as existing in discrete rotational states, with each possible state chosen to coincide with a region of low potential energy. States differ in relative energy and thus in Boltzmann-weighted probability. Discrete states are defined only around bonds that allow torsion, such as single bonds [17].
In our previous computational work on polybutadiene, we found that temperature dependence of average chain size can largely be attributed to large relative changes in low probability density for larger chains. We termed this the "taut conformation effect" [17]. This work studies the impact of such taut conformations on average sizes of polypropylene and polystyrene chains. Size properties of random chain conformations of polypropylene and polystyrene were studied by calculating the distribution of single chain conformations [17]. We assumed that multiple random conformations of a single chain in a theta solvent will have the same statistics as multiple chains generated in the bulk under the same conditions.
Thus using a single chain in each computation is justified [5] and computationally cheaper. Chain size properties were computed based on ensemble averages of single chain conformations at different chain lengths and over a range of temperatures.

Methodology
The simulations consist of generating random conformations of atactic, amorphous polypropylene and polystyrene single chains. Relative probabilities of local chain geometry guide the choice of torsion angles for each chain. Temperature ranges for polypropylene and polystyrene were selected to ensure the chains were in an amorphous state above the glass transition (T g ) [1], which is ∼258 K for atactic polypropylene [25] and ∼372 K for atactic polystyrene [26]. 100,000 chains were generated for a single chain length (degree of polymerization x = 50, 120) at temperatures T = 275, 300, 323, 343, 375, 400, and 413 K for polypropylene, and T = 375, 400, 413, 430, 450, 475, and 500 K for polystyrene. 100,000 single chains of polypropylene and polystyrene were generated at x = 15, 25, 50, 75, 100, 120, and 250 repeat units at 300 K for polypropylene and 400 K for polystyrene.
At each temperature and chain length, ensemble averages of chain size parameters were computed using an equal weighting for each chain.

Chain Geometry and Internal Energy
Details of chain start location, simulation box dimensions and total energy of each system are discussed in our previous work [17]. Bond lengths and bond angle supplements used in our computations are shown in tables 3, 4, and 5. They were obtained from Suter and Flory [11] for polypropylene, and Yoon, Sundararajan and Flory [23] for polystyrene. These are shown in tables . Table 3: Bond length geometries [11,23] Bond Length (Å) C-C 1.53 C-CH 3 1.51 C-C 6 H 5 1.51 C-H 1.09 Table 4: Bond angle geometries for polypropylene [11] Angle Table 5: Bond angle geometries for polystyrene [23] Angle The torsional angle states around C-C single bonds were based on regions of low potential energy. Suter and Flory [11] used five discrete rotational isomeric states for each torsional bond of polypropylene while Yoon, Sundararajan and Flory [23] chose two for polystyrene. The isomeric states for polypropylene correspond 11], and φ = 10 • and 110 • for polystyrene [23].

Chain Generation
Transformation matrices convert bond vectors from one reference system to another [3] and were used to determine atom positions within each single chain from internal coordinates. Two transformation matrices were used per repeat unit of polypropylene and polystyrene: one for the -CHR group and other for the -CH 2 group [11,23].
Each polymer chain was built in an atom-by-atom manner. Torsion was con-  Atoms from C i to C i+1 and their pendant atoms comprise a single repeat unit. Numbering employs Flory's convention [3]. The rotational isomeric states (i.e. torsion angles) for each repeat unit were selected randomly based on probabilities that account for relative energies of adding two consecutive torsion angles. These relative energies are statistical weights which are stored in statistical weight matrices [3]. Such weights were suggested by Suter and Flory [11] for polypropylene and Yoon, Sundararajan and Flory [23] for polystyrene. Chains were generated based on Boltzmann-weighted probabilities that account for relative energies of adding two consecutive torsion angles.
The probabilities were computed on the basis of a partition function [3]

Chain Size Parameters
Squared end-to-end distance r 2 is an important chain size parameter and was calculated as where r x , r y , r z are the x, y and z coordinates of the end-to-end distance vector r.
The vinyl polymer chains generated here are of the form CH 3 -(CHR-CH 2 ) x−1 -CHR-CH 3 where the chain start and end carbon atoms define the end-to-end distance.
Periodic boundary conditions [27] were taken into account to ensure the proper image of each atom position was employed.
The squared radius of gyration (r 2 g ) relates to the distance of each atom in the polymer chain from the center of mass, where r j = (x j , y j , z j ) T is the position vector of atom j of a polymer chain, and other terms are defined analogously. The overbar indicates an average over all atoms in a single chain conformation. The radius of gyration matrix was diagonalized and its eigenvalues were sorted in descending order (λ 1 ≥ λ 2 ≥ λ 3 ).
λ 1 corresponds to the principal direction with the longest dimension and λ 2 , λ 3 correspond to secondary directions. This approach effectively represents the size of a polymer chain using a rotational ellipsoid with a different size in each direction, rather than with a hollow sphere of radius r g that has the same mass and moment of inertia as the polymer chain [17]. The squared radius of gyration equals the sum of the three eigenvalues [28], The squared radii of gyration obtained from the eigenvalues (eq. 86) are equivalent to those from the number-weighted equation (eq. 83). Chain size parameters were quantified to analyze chain behavior at different lengths and over a range of temperatures.

Results and Discussion
A balance between intra-chain and inter-chain interactions can be inferred from the relationship between end-to-end distance and chain length. The rootmean-squared end-to-end distance increases with the number of repeat units as where ν is the scaling exponent [4]. Figure 31 shows that a scaling exponent of 0.553 best describes the simulated polypropylene data while for polystyrene it is 0.569. In both cases, ν is in between that of a theta solvent (ν = 0.5) and a good solvent (ν = 0.6) [4]. This indicates that a balance exists between intra-chain and inter-chain interactions in the simulated vinyl polymers.  Figure 31: Scaling of root-mean-squared end-to-end distance with chain length for polystyrene (blue-unfilled) at 400 K and polypropylene (red-filled) at 300 K. These symbols are used throughout unless otherwise specified. Dashed lines indicate best fits using scaling exponents of 0.569 and 0.553 for polystyrene and polypropylene, respectively.
Characteristic ratio (C n ) [3,4] can be written as where n is the number of backbone bonds along a polymer chain and l is the bond length of C-C. For vinyl polymer chains, there are two C-C bonds per repeat unit thus n = 2x. Figures 32 and 33 show characteristic ratio as a function of temperature.  Our simulation results are in decent agreement with prior model results but some disagreement is observed between our results and experimental intrinsic viscosity results. This could be the result of our chains not being in unperturbed or theta conditions as opposed to the experiments. Experimental [12,13,14] and theoretical [23,24,29] characteristic ratios were obtained for polystyrene at temperatures (∼300 K) below its glass transition temperature. The decreasing behavior seen in figs. 32 and 33 are indicative of negative temperature coefficients for both polypropylene and polystyrene.
Temperature coefficient (d ln r 2 /dT ) gives a qualitative and quantitative measure of chain size dependence on temperature. Our results may be compared to theoretical and experimental data by assuming n and l to be independent of temperature [17] d ln  [18] acknowledged that discrepancies may arise between experimental and theoretical temperature coefficients due to the sheer difficulty in obtaining theta conditions. Mark [30] mentioned about experimental difficulties with stability of swollen networks which impacts temperature coefficient calculations. It is possible that analogous difficulties could impact shrunken networks as well.       and Rapold et al. [24] used RIS models at 300 K for polystyrene and considered  [23,24] at 300 K.
The mean-squared radius of gyration r 2 g 0 should equal 1/6 of the meansquared end-to-end distance r 2 0 under unperturbed conditions and in the limit of long chains without long branches [32]. Figure 41 shows a ratio further above 6 for polystyrene than for polypropylene for the same chain length. It was also observed that r 2 / r 2 g decreased with increasing chain length but did not reach 6 indicating that the chains of 250 repeat units were not yet at the long chain limit.  Figure 41: Ratio of squared end-to-end distance to squared radius of gyration ( r 2 / r 2 g ) vs. inverse of repeat unit x for polystyrene and polypropylene chains at 400 K and 300 K respectively.

Conclusions
Rotational Isomeric State approach was used to generate random conforma-

Introduction
The single largest use of asphalt is in construction of roads and highways encompassing approximately 85% of total asphalt production [1]. Roads are paved with hot mix asphalt (HMA), which is a combination of stone, sand, or gravel held together by asphalt binder [1]. In its simplest form, asphalt binders are residues produced from crude-oil blends which are then processed with additives, usually polymers, to enhance their properties as thermoplastic adhesives for paving [2].
The lifetime of roads is closely related to the rheological properties of binders in response to deformation and temperature over time. These binders are characterized for roadway applications according to the Superpave performance-grading (PG) system from the Strategic Highway Research Program (SHRP) based on their environmental conditions of use [1]. This grading takes into account the magnitude and temperature dependence of binder viscoelastic properties [3,4].
Over time binders age, harden, and embrittle, contributing significantly to the deterioration of paved surfaces [5]. This loss of durability has been linked to shortterm aging from high temperature processing conditions as well as long-term aging during service life [6,7,8].
Rolling thin-film oven (RTFO) simulates short-term aging of binders. It is used to represent changes in binder during HMA production and placement due to oxidation [6,7,9,10]. Long-term aging is represented through use of a pressure aging vessel (PAV) to simulate in-service aging from combined effects of time, traffic and environment [6,7,9,10].
Stress-strain studies on aged binders can yield information pertaining to stiffness, resistance to fatigue and thermal cracking over different periods of time and temperature [11]. It is crucial both from a financial and practical perspective to predict asphalt durability and viscoelastic properties before investing in a material of construction meant to last a significant amount of time. Molecular modeling is a useful tool that is used to cheaply and efficiently compute and predict viscoelastic properties of asphalt binders, sometimes at time and temperature scales which are difficult and expensive to obtain experimentally.
The objective of this work was to define a self-consistent physics-based asphaltspecific model to characterize and predict mechanical behavior of binders such as those used in Rhode Island. Viscoelastic properties were obtained experimentally using a rheometer and the data was fit using time-temperature superposition [12], a technique dating back to the 1940's and 50's when Leaderman [13], Andrews [14], and Tobolsky [15] used it to predict polymer properties. Time-temperature superposition allowed predicting viscoelastic properties of binders over a broader range of time scales than experimentally measured. Modeling asphalt behavior over wide ranges of temperature, stress, and strain conditions through a well-defined physical basis allows for the prediction of macroscopic properties from microscopic behavior [16,17] with the potential of enhancing pavement designs.
Deformation measured through stress-strain behavior leads to quantifiable moduli [12,11]. As strain is repeatedly applied to a binder sample, the complex shear modulus (G * ) can be considered to be the sample's total resistance to the deformation, while phase angle δ is the lag between the applied strain and measured stress [12]. Since complex modulus (G * ) and tan δ can be related to storage (G ) and loss moduli (G ) through |G * | = G 2 + G 2 and tan δ = G /G , modeling complex modulus is equivalent to modeling storage and loss moduli [12]. Wellestablished Maxwell models [11,12] accurately predict storage and loss moduli over a continuous distribution of frequencies and relaxation times, and were used in this work. The frequency dependence in complex modulus arises from the relation between experimental time scales and binder relaxation times [12].
Relaxation time distributions can relate rheological parameters to chemical change effects on the molecular level [16]. The approach of our model is such that all viscoelastic properties can be determined in the linear viscoelastic region through a well-defined distribution of relaxation times [12,16]. Zhang and Greenfield [18,19] used this approach and identified typical molecular relaxation times in binders. Jongepier and Kuilman [20,21] proposed a log normal distribution of relaxation times. Instead of following such a Gaussian distribution, our relaxation times were restricted within finite limits [16]. Dickinson and Witt [22] computed storage and loss moduli at different frequencies and developed a relaxation spectra.
They proposed a hyperbolic relationship between complex modulus and frequency, and concluded that their spectra was not consistent with Jongepier and Kuilman's [21] assumption of log normal distribution of relaxation times.
As roads are subjected to cyclical stresses from traffic, the strain in the arrangement of binder molecules fluctuates and ultimately reaches a residual limit with time [12]. As temperatures drop, pavements contract and build up additional internal stresses which without sufficient time to relax can lead to cracking. Bending beam rheometer (BBR) measures the low temperature thermal resistance or flexibility of PAV-aged binders. In the BBR experiment, a constant load is applied to the the center of a PAV-aged binder beam for a fixed period of time, and the beam deflection and stiffness are measured as a function of time (AASHTO T 313). In this work, creep compliance computations from viscoelastic models fit to rheological data, along with stress-strain relations, allowed for prediction of low temperature thermal cracking resistance of PAV-aged binders. Additionally, molecular motion and stiffness of binders were qualitatively and quantitatively computed and compared using moduli master curves, relaxation and retardation time distributions, and stress relaxation modulus.

Experimental
The   For both strain and frequency sweep tests, the samples were allowed to equilibrate for 10 minutes prior to testing. If the test temperature was lower than 46 • C, the samples were loaded in the rheometer at 46 • C, trimmed, cooled to the test temperature, and equilibrated prior to testing. For test temperatures of 46 • C and higher, the samples were loaded at the test temperature, trimmed, equilibrated, and tested. Additionally, compliance testing was done on the PAV-aged samples.
Flexural creep strain of PAV-aged binders was determined using a Cannon Model TE-Bending Beam Rheometer (BBR) (AASHTO T 313).

Time-Temperature Superposition
Time-temperature superposition interrelates modulus, density, and temperature at times t and a T = t/t r or frequencies ω and T = ω/ω r by assuming a common functional form [12,16,23] where a T is the horizontal or frequency shift factor, T r is the reference temperature, and t r is the reference time. This technique extends beyond measurement time scales in order to predict binder mechanical properties and molecular motions at experimentally unfeasible conditions.
The same physical processes were assumed to apply across the temperature range between T and T r [12]. Density and temperature prefactors were accounted for in the time-temperature superposition by means of the vertical shift factor b T (b T = ρT /ρ r T r ) [12]. The same density was assumed across all binders and varied with temperature, leading to the same thermal expansion coefficient and b T across all binders. Combinations of a T and b T were used to shift experimentally obtained moduli, compliance and tan δ data at all temperatures to obtain master curves.
The a T were chosen based on the moduli master curves and a linear form of the Williams-Landel-Ferry (WLF) equation [24] 1 The reference temperature was chosen to be 34 • C, considering strain and frequency sweeps were done at this temperature using both 8 mm and 25 mm plates, and assigned a T and b T of 1. The same shift factors per sample were applied to generate the tan δ and compliance master curves because their time-temperature frequency shifts originate from the same mechanism [16,23].

Models for Determining Viscoelastic Properties
Maxwell models are used to represent viscoelastic properties of polymers.
These models are based on the spring and dashpot concept where spring represents the elastic segment and dashpot represents the viscous segment of the viscoelastic material. The Maxwell model assumes the spring and dashpot to be in series and is used for computing moduli (G) and stress relaxation. The Voigt model assumes the spring and dashpot to be in parallel and is used for computing creep compliance (J) [11].
Stress relaxation modulus G(t) and creep compliance J(t) from continuous and discrete models are given by [12] The plateau modulus at very low temperatures or high frequencies, where the polymer exhibits glassy behavior, is known as glassy modulus (G g ) [12]. Zeroshear viscosity (η 0 ) can be defined as the viscosity of a sample at very low shear rates when the angular frequency is almost zero [12]. Both glassy modulus and zero shear viscosity [12] were computed using the continuous Maxwell model Additionally, zero-shear viscosity can be defined by [12] η 0 = lim ω→0 dG (ω) dω (106) and determined from lowest frequency range of loss moduli master curves.

Sample Damage
Strain sweeps yield strain ranges over which binders are in the linear viscoelastic region. Additionally, they give indication of binder stability when subjected to shear. A fresh binder sample was used at each temperature over the angular frequency range. Phase angles (δ) and complex moduli (G * ) can be used to deter-mine viscoelastic properties such as storage and loss moduli (G and G ) as well as storage, loss, and complex compliance (J , J , and J * ). Stability of binders is indicated by agreement between initial and final strain sweep data at 1 and 10 rad/s respectively, with disagreement signaling sample damage. Figure 42 shows good agreement between initial and final data sets at 70 • C and disagreement at 10 • C. Over all temperatures tested, sample damage was maximum at 10 • C and no damage was observed from 34 • C upward for all samples. The linear viscoelastic region, indicated by phase angle δ being independent of strain amplitude, spans a more broad strain range at higher temperatures for all angular frequencies. Molecular motion is restricted at lower temperatures; samples are less flexible and get damaged when subjected to high shear as indicated by disagreement between initial and final data sets at 10 • C in figure 42. This is also indicated by a more pronounced change in phase angle with strain at lower temperatures [23]. Sample stability is desirable at higher temperatures considering that processing conditions such as hot mix asphalt for roads are carried out at elevated temperatures.

Time-Temperature Superposition and Viscoelastic Models
Horizontal (a T ) and vertical (b T ) shift factors were applied to experimental storage and loss moduli data relative to the reference temperature 34 • C. Details of these shift factors are given in table 11. The horizontal shift factors were modeled using the WLF equation to produce well-fit master curves. The WLF parameters for all binders are listed in table 13. Figure 43 shows decreasing slope with aging i.e. maximum slope for original or unaged binder and minimum for PAV-aged. The increasing range of a T with aging results in flatter slopes and decreased ratio of WLF parameters (C 2 /C 1 ).
The temperature dependence of horizontal shift factor a T can yield an Arrhenius relationship with 1/T through [12] where ∆H can be interpreted as the activation energy for viscoelastic relaxation and R is the universal gas constant. From table 12 and the flatter slopes in figure   43, it can be inferred that greater activation energies are required for molecular motion of aged binders.  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

% A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

% A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

% % A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Witt [16,22]. The contributions from aging processes at shorter relaxation times were similar, while greater contributions were observed for PAV-aging at longer relaxation times. This indicated that aging of binders yield additional relaxation mechanisms while retaining faster relaxations, as inferred earlier for other asphalt binders [16]. Original-2 and Original-3 showed highest contribution to modulus at lowest relaxation times and least at highest relaxation times, indicated by a sharp monotonic decrease. PAV-1 showed a slightly different shape than other binders and seemed to have peaked around ∼10 −4 seconds. Analogous to H(τ ), the contributions from aging processes at shorter relaxation times were similar while greater contributions from PAV-aging were observed at longer relaxation times for L(τ ). Using parameters from the continuous models, stress relaxation modulus G(t) and creep compliance J(t) were computed. Discrete data were in excellent agreement with continuous but were not included here. Figure 49 shows that at lower times G(t) are closely spaced for all binders and decrease with time. Original binders relax the fastest followed by RTFO-aged and then PAV-aged binders. Aging increases binder stiffness which is corroborated by PAV-aged binders having higher G(t) than RTFO-aged and original binders over the entire time scale. All binders show closely related J(t) at shorter times and increase with time due to contribution from viscous flow [12]. Original-2 and Original-3 show similar J(t) over the entire time scale while G(t) shows minor differences, with Original-2 being higher at shorter times and Original-3 higher at longer times. Low J(t) at shorter times correspond to absence of any configurational re-arrangements in the microstructure of asphalt binders [12]. Aged samples creep slower than unaged ones. G(t) can be predicted to reach zero at infinitely long times.

Bending Beam Rheometer (BBR)
J(t) was used to make a model comparison with experimental compliance data from BBR experiment. PAV-1 and PAV-2 binders were tested using BBR at -18 • C for resistance to low temperature thermal cracking. Binders were contained in beams and immersed in a cold solution bath before a constant load was applied for 240 seconds to measure beam deflection as a function of time. Calculations were averaged over two experimental runs for each binder. Maximum bending stress and strain were computed based on load, beam dimensions and deflections (AASHTO T 313). For the model comparison, compliance data were shifted to -18 • C using the WLF equation. Strain (γ(t)) was calculated where σ 0 is the maximum bending stress applied experimentally on the binders.

Conclusions
Asphalt binders from three sources were analyzed using experimental and modeling techniques. The PAV-1 binder was not polymer modified while Original-2, RTFO-2, PAV-2, and Original-3 binders were known to be modified with styrenebutadiene-styrene (SBS) co-polymer. Well-defined viscoelastic models that we employed were excellent fits to the experimental rheological data. Shift factors, horizontal or frequency (a T ) and vertical or density and temperature (b T ), were obtained by constructing viscoelastic master curves at a particular reference temperature (34 • C) to determine binder properties beyond experimentally feasible conditions. Horizontal shift factor a T varied largely with temperature due to greater variations in relaxation rates with temperature and sample preparation [16] while vertical shift factor b T showed minor variations with temperature.
WLF model and Arrhenius dependence of a T on 1/T suggested aged binders needed greater activation energies for rotational motion of their molecules. Signs of sample damage were observed across all samples at low temperatures while no such signs were observed at relatively higher temperatures. The restriction of molecular motions was greater at low temperatures thus resulting in sample deterioration with shear.
Moduli master curves indicated increase in sample stiffness with aging. This was corroborated by the relaxation time distribution and stress relaxation modulus.
Discrete and continuous models provided excellent fits to experimental results.
Glassy moduli for all samples were within typical bitumen values, except PAV-1 which was lower than expected. This could be attributed to the lack of polymer modifier within the binder. Zero-shear viscosity decreased with heating for all binders due to increased molecular motion and thus increased flexibility. Frequency ranges were not low enough to accurately predict zero-shear viscosity for PAV-2. Tan δ decreased with decreasing temperature for all binders and increasing frequency for most binders, except at certain temperatures for PAV-2 and Original-3 binders.
From relaxation and retardation time distributions, it was observed that the contributions from aging processes at shorter relaxation times were similar while greater contributions from PAV-aging were observed at longer relaxation times. At shorter times, stress relaxation moduli and creep compliance were closely spaced across all binders with original binders relaxing the fastest and PAV-aged binders relaxing the slowest. Low creep compliance at shorter times corresponded to absence of any configurational re-arrangements in the microstructure of the asphalt binders.
Compliance data from our model was shifted to -18 • C, using WLF, to compare rheological measurements for PAV-1 and PAV-2 binders with bending beam rheometer data. Our model showed good agreement with experimental results for PAV-1 over the entire time scale. For PAV-2, there was great agreement between model and experimental results at medium to longer times and slight deviation at shorter times. Both PAV-1 and PAV-2 binders showed comparable flexibilities at low temperature. Discrete rotational isomeric states were chosen corresponding to regions of potential energy minima for each torsional bond [1]. Six discrete rotational states were chosen for polybutadiene [2,3,4], five for polypropylene [5,6], and two for polystyrene [6,7,8]. Statistical weights were arranged in M ×M statistical weight matrices , where M is the number of discrete rotational isomeric states. Each repeat unit of polybutadiene has three torsional bonds while polypropylene and polystyrene have two torsional bonds, thus three statistical weight matrices were used for polybutadiene [2,3,4] and two each were used for polypropylene [5,6] and polystyrene [6,7,8] respectively.

A.2 Partition function
Partition function (z) can be defined as the sum of the unnormalized probabilities for all possible discrete rotational isomeric states [1]. where J and J are row and column matrices for chain start and end respectively, U γ , U γ , and U γ are statistical weight matrices for bonds preceding, of and fol-lowing repeat units, and x is the number of repeat units in a polymer chain [1].
Carrying out the matrix multiplication in equation A.2 provides one term in z for each possible combination of rotational isomeric states. The probability of a single conformation equals its contribution to z, divided by z.

A.3 Transformation matrices
Transformation matrices are used to transform bond vectors from one reference state to another one. According to Flory [1], the transformation matrix used to transform bond vectors from i + 1 frame to i frame (refer figures 1, 34 and 35) can be given as where θ i and φ i are bond angle supplements and torsional angles respectively.
These matrices depend on the torsional bonds in a repeat unit of a polymer chain.