Molecular simulation of nanocolloid rheology: Viscosity, viscoelasticity, and time-concentration superposition

Dinesh Sundaravadivelu Devarajan, Pouria Nourian, Gregory B. McKenna,# and Rajesh Khare* PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Molecular Simulation of Nanocolloid Rheology: Viscosity, Viscoelasticity, and Time-Concentration Superposition Department of Chemical Engineering Texas Tech University, Box 43121 Lubbock, TX 79409 − 3121 United States Corresponding author emails: greg.mckenna@ttu.edu,# rajesh.khare@ttu.edu* Tel.: (806) 834 − 4136, # (806) 834 − 0449* A particulate molecular model in which the solvent particles are considered explicitly is developed for studying the linear viscoelasticity of nanocolloidal suspensions using molecular dynamics simulations. Nanocolloidal systems of volume fractions ranging from 0.10 to 0.49 are studied. The hydrodynamics in these model systems are governed by the interparticle interactions. The volume fraction dependence of the relative zero shear viscosity exhibited by this molecular model PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. ABSTRACT is consistent with that reported in the literature experiments and simulations. Over the range of frequencies studied, the relative dynamic viscosity values follow the same qualitative trend as that seen in the literature experiments. The time-concentration superposition (TCS) principle is successfully applied to construct the viscoelastic master curves that span nine decades of frequency in the case of the elastic modulus and more than four decades of frequency in the case of the loss modulus. The TCS principle was observed to fail at high volume fractions that are near the glass transition concentration; this finding is consistent with the literature experimental and simulation observations. The volume fraction dependence of the shift factors used in the construction of the viscoelastic master curves is in good quantitative agreement with that of the viscosity of the nanocolloidal systems. Our results demonstrate that molecular simulations in conjunction with an explicit solvent model can be used to quantitatively represent the viscosity and the viscoelastic properties of nanocolloidal suspensions. Such particulate models will be useful for studying the rheology of systems whose properties are governed by specific chemical interactions. 1 Colloidal suspensions are commonly used in many applications including paints, cosmetics, pharmaceutics, and food products. The wide-spread industrial applications of these complex fluids have led to a significant interest in understanding the structure, flow, and end-use properties of such systems. In general, the colloidal suspensions, though structurally simple, exhibit a rich rheological behavior consisting of shear thinning, shear thickening, and viscoelasticity, and hence, PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. I. INTRODUCTION serve as a suitable model system for the study of rheology of complex fluids [1]. Furthermore, the colloidal suspension systems also show a complex phase behavior that is controlled by the volume fraction (𝜙) of the colloidal particles [2]. This feature of suspensions can also be used to systematically study the rheological behavior in different regions of the phase diagram by simply varying the concentration of these systems. The advancements in theoretical [3-10] and experimental [11-19] rheological techniques have had a significant impact in the understanding of the colloidal suspension rheology leading to elucidation of various rheological behaviors pertaining to shear thinning, shear thickening, viscosity, viscoelasticity, and normal stress differences. In addition, different simulation techniques are currently in practice for studying the rheological properties of colloidal suspensions. This computational effort predominantly consists of mesoscale simulation techniques that employ an implicit solvent model; examples are Brownian dynamics (BD) [20, 21] and Stokesian dynamics (SD) simulations [22, 23]. Other simulation techniques that have been used to study colloidal suspensions include the lattice Boltzmann (LB) [24], dissipative particle dynamics (DPD) [25], and stochastic rotation dynamics (SRD) [26]. These modeling techniques accurately capture the behavior of suspensions of colloidal particles in the size range of 0.2 − 1.1 𝜇𝑚. 2 the dimensions of the colloidal particles is of the order of 20 𝑛𝑚 or smaller. These suspensions are encountered in applications involving nanofillers such as titanium dioxide (TiO2), silicon dioxide (SiO2), carbon nanotubes (CNTs), and graphene, and in biological applications involving protein solutions. As the colloidal particle size becomes comparable to the solvent size, the colloid-solvent, the colloid-colloid, and the solvent-solvent intermolecular interactions play an important role in governing the behavior of these nanocolloidal suspensions [27]. Experimentally, PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. The systems of interest in this work are nanocolloidal suspensions in which at least one of a change in the structure and interparticle interactions resulting from a pH change has been demonstrated to affect the rheology of nanocolloidal suspensions of TiO2 [28] and SiO2 [29]. Similarly, the interactions between the ionic liquid (IL) molecules and the multiwalled carbon nanotubes (MWCNTs) were shown to govern the gelation behavior of the MWCNT-IL suspensions [30]. Experimental work has also indicated that the interactions between polymers – poly(vinyl alcohol) [31] or poly(ethylene glycol) [32] – and graphene oxide (GO) sheets can be used to tune the viscosity of the aqueous suspensions of GO and the polymer. The process of salt induced gelation of cellulose nanocrystal suspension is affected by the size of the cation that is added to the suspension [33]. The rheological properties of protein solutions are important for the processing operations. In a number of experimental studies, change in protein-protein interactions resulting from, for example, a change in the pH of the suspension has been shown to affect the rheological properties (e.g. viscosity, viscoelastic moduli) of the protein suspensions [34-36]. Further examples of the rheological features of the nanocolloidal systems can be found in a recent review article by Sharma et al. [37]. Given the importance of the specific chemical interactions and the local structure in determining the rheological properties of nanocolloidal systems, the solvent molecules can no longer be treated using a continuum representation in these systems. We 3 molecules explicitly in the nanocolloidal suspensions. Molecular dynamics (MD) simulation technique in conjunction with the explicit particulate solvent model has the unique ability to account for these chemical interactions and the structural heterogeneities in the system. MD simulations have been used to elucidate the mechanism of aqueous phase stabilization PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. thus follow the approach of using a particulate nanoscale molecular model to treat the solvent of carbon nanotube (CNT) dispersions [38] and boron nitride nanosheet (BNNS) dispersions [39, 40]. Previous MD simulation studies of colloid rheology have focused on investigating the effect of different nanoparticle shapes on the rheological behavior of suspension systems [41] as well as to obtain the shear viscosity of such systems [42]. However, the linear viscoelastic behavior of nanocolloidal suspensions has not yet been extensively studied using MD simulations. Furthermore, like experimental rheology, simulations are restricted to a limited range of frequency or time scale window. The goals of this work are two-fold: First goal is to characterize both the steady shear rheology and the oscillatory shear rheology of nanocolloidal model systems over a wide range of volume fractions and frequencies. Checking the applicability of the time- concentration superposition (TCS) principle for creating the master curves of viscoelastic moduli values and thus extending the range of frequencies accessible to simulations is another major goal of this work. The rheological properties of the model suspension systems determined in this work were seen to be consistent with the literature experimental, simulation, and theoretical studies for the range of volume fractions studied. Furthermore, the master curves of the elastic (𝐺 ′ ) and the loss (𝐺 ′′ ) modulus values constructed using the TCS principle resulted in an extended dynamic rheological map that spans nine decades of frequency in the case of 𝐺 ′ and more than four decades of frequency in the case of 𝐺 ′′ . The TCS principle was found to work only in the 𝛼-relaxation region and the loss modulus values do not superpose in the 𝛽-relaxation region. Our work 4 nanocolloidal suspensions using a particulate model for the solvent. The rest of the paper is organized as follows: In Sec. II, we discuss the details of the molecular model of nanocolloidal suspension systems and the simulation technique used to study these models. In Sec. III, the results and discussions of the relative zero shear viscosity (𝜂𝑟 ), the PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. demonstrates the ability of MD simulations for capturing the linear viscoelastic behavior of elastic (𝐺 ′ ) and the loss (𝐺 ′′ ) modulus, and the relative dynamic viscosity (𝜂𝑟∗ ) of nanocolloidal suspensions for different volume fractions (𝜙) are presented. Further, the results of the master curves of the linear viscoelastic modulus and the physical significance of the shift factors are discussed. We close our paper with a brief summary of the findings in Sec. IV. II. MOLECULAR MODEL AND SIMULATION DETAILS The simulation systems of interest in this study consisted of colloidal particles that were dispersed in an explicit solvent system. A 5: 1 colloid-solvent size ratio was employed and suspensions with colloidal volume fractions (𝜙) ranging from 0.10 to 0.49 (i.e., covering the regimes of dilute, moderate, and concentrated suspensions) were studied. A particulate model of the colloidal suspension system was used in which two components are present: the solvent particles and the colloidal particles. The background solvent particles were modeled as beads and the colloidal particles were treated as large smooth spheres. The interactions between all of the beads were modeled by the purely repulsive Lennard-Jones (LJ) potential, known as the Weeks-ChandlerAndersen (WCA) potential [43]. In particular, the solvent particles interacted with each other via the WCA potential as follows: 𝜎 12 4𝜀 [( ) 𝑟 𝑈𝑊𝐶𝐴 (𝑟) = { 𝜎 6 1 − ( ) + ] , 𝑟 ≤ 𝑟𝑐 = 21/6 𝜎 𝑟 5 4 0, 𝑟 > 𝑟𝑐 (1) 𝑡 in this work are reported in the reduced LJ units defined by 𝑡 ∗ = (𝑚𝜎2⁄ 𝑟 𝜎 𝜀 )1⁄2 , 𝜌∗ = 𝜌𝜎 3 , and 𝑟 ∗ = , etc. The quantities in the reduced units can be converted to those in the real units by expressing them in terms of the parameters of the coarse-grained particulate model used here. For example, if each of the solvent particles is mapped to a single water molecule using the TIP4P/2005 model PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. where 𝜀 and 𝜎 are the LJ parameters, and 𝑟𝑐 is the cut-off distance. All of the quantities discussed [44], then a reduced shear rate of 𝛾̇ ∗ = (𝜀⁄ 𝜂 (𝑚𝜀⁄𝜎 4 )1⁄2 =1 respectively. will correspond to 𝛾̇ 𝑚𝜎 2 )1⁄2 values of = 0.05 and a reduced viscosity of 𝜂 ∗ = 3.28 × 1010 𝑠 −1 and 0.062 × 10−3 𝑃𝑎. 𝑠, Furthermore, when mapped to a different coarse-grained model such as the MARTINI model [45] in which each of the solvent particles corresponds to four water molecules, 𝛾̇ ∗ = 0.05 and 𝜂 ∗ = 1 will correspond to values of 2.80 × 1010 𝑠 −1 and 0.143 × 10−3 𝑃𝑎. 𝑠, respectively. The timescale accessible to simulation will increase with the increasing level of coarse-graining. For simplicity, the reduced units will be represented without the asterisk symbol in the rest of this paper. The interactions between the colloidal beads and between the colloid and the solvent beads were modeled using the distance-shifted WCA potential [46, 47] as follows: 𝜎 12 4𝜀 [( ) 𝑟−𝑎 𝑈𝑠ℎ𝑖𝑓𝑡𝑒𝑑−𝑊𝐶𝐴 (𝑟) = { −( 𝜎 𝑟−𝑎 6 1 ) + ] , 𝑟 − 𝑎 ≤ 𝑟𝑐 = 21/6 𝜎 4 0, 𝑟 − 𝑎 > 𝑟𝑐 (2) In both of the cases discussed above, values of 𝜀 = 1 and 𝜎 = 1 were used. The parameter 𝑎 takes into account the size of the colloidal particles. In our simulations, the value of 𝑎𝑐𝑐 = 𝑐 𝑐 2𝑅𝑏𝑎𝑟𝑒 − 𝜎 = 4 was used for the colloid-colloid interactions, where 𝑅𝑏𝑎𝑟𝑒 = 2.5 represents the bare radius of the colloidal particles. In order to account for the depletion forces that can lead to the aggregation of colloidal particles [24], the parameter 𝑎𝑐𝑠 for the colloid-solvent interactions 𝑐 𝑠 was set to a value that is slightly smaller than its expected value (𝑎𝑐𝑠 = 𝑅𝑏𝑎𝑟𝑒 + 𝑅𝑏𝑎𝑟𝑒 − 𝜎 = 2), 6 PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. FIG. 𝟏. Snapshots of the sample simulation boxes displaying the molecular models of the nanocolloidal dispersions for volume fractions, 𝜙 = 0.20, 0.35, and 0.49. The colloidal particles that are modeled as large smooth spheres are shown in cyan color and the solvent particles (greatly reduced in size for the sake of clarity) are shown in gray color. FIG. 𝟐. Struik’s protocol for carrying out the physical aging simulations using oscillatory perturbations (Figure adapted with permission from Ref. [52]. Copyrighted by the American Physical Society 2016). 7 𝜎 2 represents the bare radius of the solvent particles; this modification allows the solvent particles to penetrate in the gap between any two colloidal particles that are approaching 𝑐 𝑠 each other at short distances. Thus, in this work, the value of 𝑎𝑐𝑠 = 0.9 × (𝑅𝑏𝑎𝑟𝑒 + 𝑅𝑏𝑎𝑟𝑒 − 𝜎) = 1.8 was used for the colloid-solvent interactions [48]. A periodic cubic simulation box of edge length 𝐿 = 50 was used for all of the molecular PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 𝑠 where 𝑅𝑏𝑎𝑟𝑒 = dynamics (MD) simulations carried out in this work. The volume fraction of the colloidal suspension systems in the central simulation box was defined as 𝜙 = 𝑁𝑐 4 𝑉 3 𝑐 𝜋(𝑅𝑏𝑎𝑟𝑒 )3 , where 𝑁𝑐 is the number of colloidal particles and 𝑉 is the simulation box volume. The number of colloidal particles varied from 𝑁𝑐 = 191 to 936 for volume fractions, 𝜙 = 0.10 − 0.49, that are studied in this work. As the colloidal volume fraction was increased, the proportionate number of solvent particles were removed from the system to maintain the mass density constant at 𝜌 = 0.60 for all of the volume fractions studied. The mass of the colloidal particles was calculated as 𝑚 = 4 𝑐 𝜌 × 𝜋(𝑅𝑏𝑎𝑟𝑒 )3 = 39.27 to ensure that the density of the colloids is the same as that of the 3 background solvent medium. Snapshots of the sample simulation boxes displaying the particulate molecular model of the nanocolloidal suspension system for volume fractions, 𝜙 = 0.20, 0.35, and 0.49, are shown in Fig. 1. The linear rheological properties of the model colloidal systems were calculated using methodology similar to that employed in our recent simulation study [49] of the rheology of atomistically detailed models of asphalt and involved applying the non-equilibrium MD (NEMD) simulation technique [50]. Five independent replicas were used to obtain all of the rheological quantities presented in this work. The particulate model systems were subjected to a steady shear strain at different shear rates (𝛾̇ ) using the SLLOD equations of motion and the sliding-brick 8 viscosity (𝜂) was obtained. A time-step of Δ𝑡 = 0.003 was used for the steady shear simulations. The systems were equilibrated for a period of 40 × 106 time-steps and this stage was then followed by a production run of 60 × 106 time-steps. It was verified that beyond the equilibration stage, i.e. in the production stage, the shear viscosity values did not change with time for all volume fractions (𝜙) studied. In the case of oscillatory shear simulations, the only difference is the use of an oscillatory shear strain that was imparted on the model systems. The resulting in-phase and PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. periodic boundary conditions [50], and the shear stress was measured from which the shear out-of-phase components of the oscillatory shear stress were then used to determine the elastic (𝐺 ′ ) and the loss (𝐺 ′′ ) modulus values of the colloidal systems. For these simulations, a time-step of Δ𝑡 = 0.002992 − 0.003021 (a different time-step value was used at each frequency (𝜔) to ensure that the oscillation time period is an integer multiple of the time-step) was used. Further, the model systems were initially allowed to equilibrate by subjecting them to an oscillatory perturbation for a period of 50 oscillations at all frequencies in the chosen range of 𝜔 = 0.004 − 0.25. Following this, in the production stage, the oscillatory shear stress response was collected for a duration of 100 oscillations in the case of lower frequencies (𝜔 < 0.015) and for a duration of 150 oscillations in the case of higher frequencies (𝜔 ≥ 0.015). Large number of oscillations were required in the equilibration and the production stages to ensure that the uncertainties in the modulus values were small. The long relaxation times of the system at high volume fractions (𝜙 ≥ 0.46) necessitated special considerations in the determination of the 𝐺 ′ and the 𝐺 ′′ values at these 𝜙. To further investigate this aspect, the Struik’s protocol [51, 52] that is commonly used to account for physical aging in creep experiments was followed at the highest volume fraction of 𝜙 = 0.49. A schematic of the modified Struik’s protocol that was employed to account for the physical aging at this 𝜙 9 the aging time or the waiting time (𝑡𝑤 ) of the model systems was at least ten times the time (𝑡𝑖 ) for which the systems are perturbed during the oscillatory shear simulations. Three sequences of Struik’s protocol were performed in the frequency range of 𝜔 = 0.04 − 0.25. We note that these aging simulations required to complete the Struik’s protocol were computationally expensive and accounted for a total duration of more than 109 time-steps. For all of the simulations reported in PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. using oscillatory perturbations is shown in Fig. 2. The Struik’s protocol was followed such that this work, the system temperature was maintained at 𝑇 = 1 using the Nosé-Hoover thermostat [53]. The MD simulations at each volume fraction were performed at the conditions of constant number of particles, volume, and temperature (constant NVT) using the LAMMPS package [54]. III. RESULTS AND DISCUSSIONS A. Molecular structure in nanocolloidal suspensions For the purpose of quantifying the structure present within the colloidal suspensions, the radial distribution functions (RDFs) of the colloid-colloid and the colloid-solvent pairs were calculated. The RDF (𝑔(𝑟)) quantifies the probability of finding two particles at a given distance normalized by the probability of occurrence of those particles at that distance in the case of a uniform distribution of the particles. Figs. 3(𝑎) and 3(𝑏) show the colloid-solvent and the colloid-colloid RDFs for volume fractions, 𝜙 = 0.15, 0.35, 0.45, 0.46, and 0.49. As seen from the figures, the magnitude of the first peak increases with an increase in the volume fraction indicating higher tendency for the short-range packing of particles at large 𝜙 values. Furthermore, 𝑔(𝑟) attains a value of 1 at a shorter distance in the case of the dilute suspension (𝜙 = 0.15) compared to that for the concentrated suspensions (𝜙 = 0.45, 0.46, and 0.49), where the 𝑔(𝑟) exhibits peaks at distances of the order of 4 particle diameters. This long-range order in the structure indicates that 10 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 (𝒂) (𝒃) FIG. 𝟑. Radial distribution functions (𝑔(𝑟)) of (𝑎) the colloid-solvent and (𝑏) the colloid-colloid pairs for volume fractions, 𝜙 = 0.15, 0.35, 0.45, 0.46, and 0.49. 11 system at these high volume fractions. Similar observations were also reported in the theoretical studies of hard sphere suspensions [55]. As seen from Fig. 3(𝑎), the colloid-solvent RDF shows a strong first peak at a distance of 2.85 at all volume fractions. We note that this peak position is approximately the same as 𝑎𝑐𝑠 + 𝜎 (= 1.8 + 1), indicating that the first solvent layer is in close contact with the colloidal particles. The RDF plot also shows a smaller second peak at a distance PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. the colloidal systems are undergoing a transition from a dilute fluid system to a densely-packed of 4.05 (≈ 𝑎𝑐𝑠 + 2𝜎) which is indicative of a weak tendency for the formation of a second layer of solvent beads around the colloidal particles. Fig. 3(𝑏) shows a first peak in the colloid-colloid RDF at a distance of 5.5 (≈ 2 × [first peak location in colloid-solvent RDF]). Thus, this peak location suggests that two colloidal particles are not in direct contact but are separated by a layer of solvent atoms, as was also suggested by the first peak in Fig. 3(a). Furthermore, in connection with the observation of a small second peak in the colloid-solvent RDF, there is also a shoulder in the colloid-colloid RDF. This observation is consistent with the configurations where colloidal particles are separated by two layers of solvent particles. We note that the observed solvent layering around the colloidal particles results in a larger separation between them, which in turn, affects the suspension structure and hence, the rheology of these systems, as compared to that observed in an implicit solvent model system [56]. A weak second peak in the colloid-colloid RDF is only observed for volume fractions 𝜙 ≥ 0.35, indicating that the colloids exhibit only very short range packing at the lower volume fractions. Furthermore, the crystallization effects were not observed at these 𝜙 values as indicated by the absence of peaks as those expected, for example, when an FCC crystal is formed in colloidal systems [57]. Following the work by Kohale and Khare [58], the hydrodynamic radius of the colloidal particles is deduced from the location of the first peak in the colloid-solvent RDF (see Fig. 3(𝑎)) and is found to have a value of 𝑅𝐻𝑐 = 2.85 12 B. Steady shear viscosity As a first step towards studying the rheology of colloidal dispersions using molecular simulations, the rheological properties of the particulate molecular model that is used in this study were compared with the prior literature reports of the rheological behavior of the colloidal suspension systems [59]. For this purpose, the NEMD simulation technique was used to determine the relative PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. at all volume fractions. zero shear viscosity (𝜂𝑟 ) values as a function of the colloidal volume fraction (𝜙). The model systems were subjected to steady shear at different shear rates (𝛾̇ ) and the corresponding shear stress component (𝜏𝑥𝑦 ) was measured; the shear viscosity values (𝜂) as a function of the shear rate (𝛾̇ ) were then determined as 𝜂 = − 𝜏𝑥𝑦 𝛾̇ (see Fig. 4). As expected, the viscosity values increase with an increase in the volume fraction of the suspension. In order to test the effect of the interparticle interactions on the viscosity values, for one of the volume fractions (𝜙 = 0.40), an additional set of simulations was carried out by changing the interaction between only the colloid and the solvent particles from WCA to full LJ potential with the value of well-depth set to be 𝜖 = 2. As seen, this simple change in the colloid-solvent interaction resulted in an approximately factor of two increase in the viscosity of the suspension compared to the value obtained using the purely repulsive LJ or WCA interactions between the colloid and the solvent particles. The colloidal suspensions show Newtonian behavior at low volume fractions (𝜙 ≤ 0.45) while a strong shear thinning behavior is observed at the high volume fractions (𝜙 ≥ 0.46) that were studied in this work. The value of zero shear viscosity (𝜂0 ) was determined from the Newtonian plateau observed at the low shear rates. In the case of high volume fractions (𝜙 ≥ 0.46), the Newtonian plateau could not be attained for the shear rates studied (down to 𝛾̇ = 13 PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. FIG. 𝟒. Shear viscosity (𝜂) as a function of the shear rate (𝛾̇ ) for colloidal suspensions of different volume fractions (𝜙). We note that the hollow green squares are the viscosity results for 𝜙 = 0.40 that were obtained by modeling the colloid-solvent interactions via the LJ potential (LJcs ) with the value of well-depth set to be 𝜀 = 2. 14 plateau at this 𝜙 are too expensive computationally. For example, the simulations to obtain the 𝜂 value at the lowest 𝛾̇ studied here, 𝛾̇ = 1 × 10−5 took about 3 days and we thus expect that the simulations at a 𝛾̇ that is an order of magnitude lower will take about 3 − 4 weeks to attain a reliable 𝜂 value. The relative zero shear viscosity values at different volume fractions were then PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 1 × 10−5 ). Simulations at even lower shear rates that are required to capture the Newtonian calculated as 𝜂𝑟 = 𝜂0 𝜂𝑠 , where 𝜂𝑠 represents the zero shear viscosity of the pure solvent system. The zero shear viscosity of the solvent (at 𝜌 = 0.60) as calculated using the NEMD simulations in this work is 𝜂𝑠 = 0.819 ± 0.122. The 𝜂𝑟 values of the colloidal suspensions so-obtained are then compared with those reported in the literature experimental [60-63] and simulation [64-66] studies of hard-sphere suspensions. As can be seen from Fig. 5(𝑎), after an initial slow increase in 𝜂𝑟 with an increase in the volume fraction, the simulated 𝜂𝑟 values rise rapidly at high volume fractions (𝜙 > 0.40). We note that the model colloidal particles used in this study interact via a purely repulsive LJ potential and thus our model systems are not exactly the same as the hard sphere suspensions. Despite this difference, it is seen that the 𝜂𝑟 values obtained from the NEMD simulations are in good agreement with the literature experimental and simulation results. The only exception is that observed at a high volume fraction value of 𝜙 = 0.45, for which the 𝜂𝑟 value obtained in this work is slightly higher than the literature value. However, our viscosity result at 𝜙 = 0.45 is in good agreement with that reported from the MD simulations by in’t Veld et al. [42] which employed a different interaction potential for the colloid-solvent and the colloid-colloid interactions. Furthermore, it is also seen that at low volume fractions (𝜙 ≤ 0.20), the simulated 𝜂𝑟 values follow the Einstein viscosity equation, 𝜂𝑟 = 1 + 2.5𝜙. Thus, our model colloidal suspension systems with colloid to solvent size ratio of 5: 1 exhibit behavior that is consistent with the continuum expectations. 15 PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. (𝒂) (𝒃) FIG. 𝟓. (𝑎) Comparison of the relative zero shear viscosities (𝜂𝑟 ) of colloidal suspensions as obtained from the NEMD simulations (shown as red diamonds) with those that are reported in the literature experimental and simulation studies. The dashed orange line and the solid dark magenta line are the predictions obtained by fitting the Krieger-Dougherty (KD) equation and the model based on the mode coupling theory (MCT) to the 𝜂𝑟 values obtained from NEMD. Best fit values of 𝜙𝑚 = 0.476 and 𝜙𝑚 = 0.519 were obtained in the case of KD equation and MCT model, respectively. The dash-dot black line is the prediction of the 𝜂𝑟 values using the Einstein viscosity equation for suspensions. (𝑏) same as (𝑎) but all of the 𝜂𝑟 values from the literature studies and the simulated NEMD values of 𝜂𝑟 are rescaled using the 𝜙𝑔,𝑀𝐶𝑇 value obtained from the respective MCT fit to each set of 𝜂𝑟 values. 16 by semi-empirical models. The two commonly used models are the Krieger-Dougherty (KD) equation [67] and the model based on the mode-coupling theory (MCT) [59]. The KriegerDougherty equation can be written as, 𝜂𝑟 = [1 − PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. The volume fraction dependence of viscosity of the colloidal systems is often represented 𝜙 −2.5𝜙𝑚 𝜙𝑚 ] (3) where 𝜙𝑚 is the maximum packing fraction of the colloidal particles and [𝜂] = 2.5 is the intrinsic viscosity of the colloidal particles. As seen in Fig. 5(𝑎), a value of 𝜙𝑚 = 0.476 is found to fit the 𝜂𝑟 values from the NEMD simulations. The semi-empirical model based on the mode coupling theory (MCT) can be written as ′ 𝜂−𝜂∞ 𝜂𝑠 = 0.9𝜙 2 [1 − 𝜙 −2.46 𝜙𝑚 ] (4) ′ where 𝜂∞ represents the high frequency viscosity that is obtained by using a semi-empirical correlation as [68], ′ 𝜂∞ 𝜂𝑠 = 1+1.5𝜙(1+𝜙−0.189𝜙2 ) 1−𝜙(1+𝜙−0.189𝜙2 ) (5) Using Eqs. (4) and (5), it was shown that in the case of hard sphere suspensions, the divergence of viscosity occurs at the glass transition value given by 𝜙𝑚 = 𝜙𝑔 = 0.57 [59]. However, when the MCT based model was used to fit the 𝜂𝑟 values obtained from the NEMD simulations in this work, a value of 𝜙𝑚 = 0.519 is obtained. Thus, the fits of the simulation results to the Krieger- Dougherty and the MCT based models indicate that the glass transition concentration (𝜙𝑔 ) for our particulate colloidal model systems is in the range of 0.476 − 0.519. We note that, on the theoretical side, Bengtzelius et al. [69] have reported a value of 𝜙𝑔 = 0.516 in the case of a hard 17 using the 𝜙𝑔,𝑀𝐶𝑇 value obtained by fitting the MCT model to each set of 𝜂𝑟 values individually (see Table I) and a plot of 𝜂𝑟 vs. rescaled 𝜙 (= 𝜙 𝜙𝑔,𝑀𝐶𝑇 ) is shown in Fig. 5(𝑏). As seen in the figure, the average 𝜂𝑟 values obtained from the NEMD simulations are slightly lower than the literature 𝜂𝑟 values at low 𝜙 values. This difference can be attributed to the use of a soft potential PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. sphere model system. We have also rescaled the 𝜙 values reported in this work and in the literature for simulating our model systems that leads to a steep increase in the 𝜂𝑟 value at 𝜙 = 0.45 (see Fig. 5(𝑎)). TABLE 𝐈. Values of 𝜙𝑔,𝑀𝐶𝑇 that were obtained by fitting the MCT model to the relative zero shear viscosity (𝜂𝑟 ) values reported in this work and in the literature studies. NEMD simulations (this work) 𝝓𝒈,𝑴𝑪𝑻 in’ t Veld et al. (MD simulation) [42] 0.524 Sriram et al. (experiment) [60] 0.571 Segre et al. (experiment) [61] 0.553 Cheng et al. (experiment) [62] 0.568 Shikata and Pearson (experiment) [63] 0.616 Foss and Brady (SD simulation) [64] 0.571 Phung (SD simulation) [65] 0.574 Pan et al. (DPD simulation) [66] 0.652 Source 18 0.519 In the previous section, it was shown that for our model systems up to 𝜙 = 0.45, the volume fraction dependence of viscosity is quantitatively similar to that observed in the experiments and in the previous simulation studies. In this section, the results of the elastic (𝐺 ′ (𝜔)) and the loss (𝐺 ′′ (𝜔)) modulus values of colloidal systems for volume fractions ranging from 𝜙 = 0.15 − 0.49 as obtained using the oscillatory NEMD simulation technique are presented. The shear stress PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. C. Linear viscoelasticity response of the systems over different oscillation periods was averaged and then fitted using the relationship, 𝜎𝑥𝑦 = 𝜎0 sin(𝜔𝑡 + 𝛿), to determine the shear stress amplitude (𝜎0 ) and the phase lag (𝛿). The linear viscoelastic modulus (𝐺 ′ and 𝐺 ′′ ) values were then calculated from the in-phase and the out-of-phase components of the stress response as 𝐺 ′ = 𝜎0 𝛾0 cos(𝛿) and 𝐺 ′′ = 𝜎0 𝛾0 sin(𝛿). For the purpose of ensuring that the linear viscoelastic modulus values are not affected by the amplitude of the shear strain (𝛾0 ) imparted on the simulation box, an initial set of oscillatory simulations was performed to obtain the 𝐺′ and the 𝐺 ′′ values as a function of 𝛾0 at the frequency of 𝜔 = 0.1. Fig. 6 shows the strain amplitude (𝛾0 ) dependence of the elastic (𝐺 ′ ) and the loss (𝐺 ′′ ) modulus values that were obtained at a frequency of 𝜔 = 0.1 for different colloidal volume fractions (𝜙). As can be seen from the figure, 𝐺 ′′ values exhibit a plateau region (indicative of linear response) at low 𝛾0 values, while at higher values of strain amplitudes, both 𝐺 ′ and 𝐺 ′′ values show a dependence on 𝛾0 , thus indicating that the system response is non-linear in this region. In order to reduce the uncertainties in the results, the simulations for determining modulus values as a function of frequency at each 𝜙 were carried out at the highest possible 𝛾0 value that resides in the linear region. We note that since the elastic modulus (𝐺 ′ ) values are of very small magnitude for 𝜙 ≤ 0.30, a clear plateau region was not observed for these, and hence, only 𝐺 ′′ 19 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 (𝒂) (𝒃) FIG. 𝟔. (𝑎) Elastic (𝐺 ′ ) and loss (𝐺 ′′ ) modulus values as a function of the strain amplitude (𝛾0 ) at a a frequency of 𝜔 = 0.1 for colloidal suspensions of different volume fractions (𝜙). 20 highest 𝛾0 value at which the linear response is observed at each 𝜙 decreases with an increase in the volume fraction; similar observation was also made in an experimental study by Mason and Weitz [11]. The linearity of the system response to the oscillatory shear was further tested by taking a Fourier transform of the stress signal at each frequency and verifying that only a single peak, that at the input frequency, was observed. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. values were used to select the 𝛾0 value to be used in simulations for modulus determination. The For the system at the highest volume fraction (𝜙 = 0.49), the Newtonian plateau required for obtaining the zero shear viscosity (𝜂0 ) was not attained at the shear rates studied in our work, thus indicating long relaxation times at 𝜙 = 0.49. Given this observation, Struik’s protocol as shown in Fig. 2 was followed to check for any possible physical aging effects at 𝜙 = 0.49 (see Fig. 7). For each of the five replicas, the modulus values at the lowest frequency of 𝜔 = 0.04 were first obtained and the resulting configuration from these oscillatory shear simulations was then used to obtain the modulus values for these replicas at the next lowest frequency and so on. It is observed from the figure that the modulus values are the same for the three sequences of Struik’s protocol. We note that if such a procedure of determining the modulus values in a sequential fashion from the lowest to the highest frequency is not followed (i.e., if the same starting structure is used to obtain the modulus values at different frequencies in parallel), the modulus values obtained in different sequences of Struik’s protocol are different. We believe that this observation is an indication that the system undergoes physical aging at this 𝜙 but quickly reaches its equilibrium aged state by the end of the low frequency simulation run. Based on this observation, for 𝜙 = 0.455 − 0.49, the elastic (𝐺 ′ ) and the loss (𝐺 ′′ ) modulus values were obtained from the sequential MD simulation runs (i.e., lowest to highest frequency) over the frequency range of 𝜔 = 0.004 − 0.25. 21 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 (𝒂) (𝒃) FIG. 𝟕. (𝑎) Elastic (𝐺 ′ (𝜔)) and (𝑏) loss (𝐺 ′′ (𝜔)) modulus values before and after the application of Struik’s protocol for 𝜙 = 0.49. 22 Having established that the system response is linear in our simulations and that the modulus values at the highest volume fraction are free from the physical aging effects, we now focus on the frequency and the volume fraction dependence of the viscoelastic modulus values. The elastic (𝐺 ′ ) and the loss (𝐺 ′′ ) modulus values as obtained from the NEMD simulations for volume PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 1. Frequency dependence of viscoelastic modulus fractions, 𝜙 = 0.15 − 0.49 are shown in Figs. 8(𝑎) and 8(𝑏), respectively. As expected, both the 𝐺 ′ and the 𝐺 ′′ values of the colloidal suspensions increase with an increase in the frequency and an increase in the volume fraction up to 𝜙 = 0.45. At high volume fractions of 𝜙 ≥ 0.455, a plateau in 𝐺 ′ is observed whereas at these 𝜙 values, the 𝐺 ′′ values continue to monotonically increase as a function of the frequency (𝜔). Furthermore, the value of the ratio 𝐺′ 𝐺 ′′ increases with an increase in the 𝜙 value such that at 𝜙 = 0.45, the ratio has a value close to 1 at all frequencies studied while the ratio becomes greater than 1 (and shows a decrease in the frequency) for 𝜙 ≥ 0.455 (see Fig. 8(c)). This observation of an increase in the elastic behavior accompanied by a plateau in 𝐺 ′ values at high volume fractions (𝜙 ≥ 0.455 ) signals the approach to the glass transition concentration (𝜙𝑔 ) [11]. These findings are also consistent with the deduction of the glass transition concentration range of 0.476 < 𝜙𝑔 < 0.519, as obtained from the semi-empirical model fits to the viscosity values. 2. Dynamic viscosity: Comparison with experiments and validity of Cox-Merz rule The dynamic viscosity values (𝜂𝑟∗ ) obtained in this work from the NEMD simulations for different volume fractions (𝜙) were compared with those obtained from the literature oscillatory active 23 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 (𝒂) (𝒄) FIG. 𝟖. (𝑎) Elastic (𝐺 ′ (𝜔)) modulus, (𝑏) loss (𝐺 ′′ (𝜔)) modulus, and (𝑐) suspensions as obtained from the NEMD simulations for different volume fractions (𝜙). 24 𝐺′ 𝐺 ′′ (𝒃) ratio of the colloidal frequency or the oscillation rate (𝛼) defined as, 𝛼= 𝑐 𝜔(𝑅𝑏𝑎𝑟𝑒 )2 2𝜋𝐷𝑐 = 𝑐 ((𝑅𝑏𝑎𝑟𝑒 )2 ⁄𝐷𝑐 ) (6) (2𝜋⁄𝜔) where 𝐷𝑐 is the diffusivity of the colloidal particles. It can be seen from Eq. (6) that the oscillation 𝑐 rate is the ratio of timescale of the diffusion of colloidal particles ((𝑅𝑏𝑎𝑟𝑒 )2 ⁄𝐷𝑐 ) to the timescale PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. microrheology experiments [60]. The comparison was carried out using the dimensionless of the oscillatory deformation (2𝜋⁄𝜔). We note that, in our work, all of the reported material parameters and 𝜔 are already in the dimensionless LJ units. The dynamic viscosity values were ′ ′′ obtained from the 𝐺 and the 𝐺 values by using the relation 𝜂 ′′ = 𝐺′ 𝜔 𝜂𝑟∗ = √𝜂′ 2 +𝜂′′ 2 𝜂𝑠 , where 𝜂 ′ = 𝐺 ′′ 𝜔 and . Fig. 9(𝑎) shows a comparison of the 𝜂𝑟∗ values obtained from the oscillatory NEMD simulations and the oscillatory active microrheology experiments [60] at different volume fractions. We note that in ref. [60], the 𝜂𝑟∗ values were presented by rescaling the volume fraction values using the effective radius of the colloidal particles; we have compared our results with the results that were reported in ref. [60] using these rescaled volume fraction values. Though the 𝜙 values used in the experiments and in our simulations are slightly different, the 𝜂𝑟∗ values obtained from the two techniques are semi-quantitatively consistent and both exhibit a Newtonian plateau at low frequencies followed by the onset of frequency thinning with an increase in the value of 𝛼. Furthermore, as expected, the dynamic viscosity (𝜂𝑟∗ ) values increase with an increase in the volume fraction (𝜙). It has been previously reported in the literature that the Cox-Merz rule completely breaks down [70, 71] or fails only at high volume fractions [60] in the case of the colloidal suspension systems. For the purpose of comparing the results of the steady shear and the oscillatory shear 25 PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. (𝒂) (𝒃) FIG. 𝟗. (𝑎) Comparison of the relative dynamic viscosity (𝜂𝑟∗ ) values obtained from the oscillatory NEMD simulations and the oscillatory active rheology experiments [60] for different volume fractions (𝜙) and (𝑏) comparison of the relative steady shear viscosity (𝜂𝑟 ) and the relative dynamic viscosity (𝜂𝑟∗ ) values to check for the validity of the Cox-Merz rule. 26 𝜂𝑟 (𝛾̇ ) = 𝜂𝑟∗ (𝜔⁄2𝜋) (7) In order to arrive at a meaningful comparison, the shear rate (𝛾̇ ) was also represented in the form of a Péclet number as 𝑃𝑒 = 𝑐 𝛾̇ (𝑅𝑏𝑎𝑟𝑒 )2 𝐷𝑐 . Fig. 9(𝑏) shows the comparison of the relative steady shear (𝜂𝑟 ) and the relative dynamic viscosity (𝜂𝑟∗ ) values obtained from simulations for our model systems PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. simulations obtained in this work, the Cox-Merz relationship was employed as follows [72]: used in this study. As can be seen from the figure, the 𝜂𝑟∗ values are quantitatively consistent with the shear viscosity (𝜂𝑟 ) values only up to 𝜙 = 0.40 for our model systems and the agreement between the two deteriorates for 𝜙 ≥ 0.45. D. Time-concentration superposition principle: Master curves of 𝑮′ and 𝑮′′ Theoretical work by Brady [73] has indicated that the real and the imaginary parts of the complex viscosity of suspensions of different volume fractions can be superposed on a master curve by scaling the frequency using the concentration dependent time scale for diffusion. Furthermore, it has also been shown that the experimental modulus values for colloidal glasses at different volume fractions (𝜙) can be collapsed onto a master curve using the time-concentration superposition (TCS) principle [74, 75]. However recently, Peng et al. [76] have argued that the TCS principle works for colloidal suspensions only when the 𝛼-relaxation and the 𝛽-relaxation processes are considered individually but fails when applied to a combination of both of these mechanisms. In order to test the applicability of this superposition principle to our model systems, we followed a similar approach by applying the TCS principle to the simulated values of the elastic (𝐺 ′ ) and the loss (𝐺 ′′ ) modulus to construct the master curves of the modulus values. We note that the solvent contribution (𝜔𝜂𝑠 ) is not removed from the 𝐺 ′′ values reported in this work. This is because the time-concentration superposition (TCS) principle is applicable only to the 𝐺 ′ and the 𝐺 ′′ values 27 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 FIG. 𝟏𝟎. Master curve of the loss tangent (tan 𝛿 ) as a function of the scaled frequency (𝑎𝜙 𝜔). The chosen reference volume fraction is 𝜙 = 0.40. 28 PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. FIG. 𝟏𝟏. Master curves of the scaled elastic (𝑏𝜙 𝐺 ′ ) and the scaled loss (𝑏𝜙 𝐺 ′′ ) modulus as a function of the scaled frequency (𝑎𝜙 𝜔). The chosen reference volume fraction is 𝜙 = 0.40. FIG. 𝟏𝟐. Van Gurp-Palmen plot of the scaled complex modulus (𝑏𝜙 ȁ𝐺 ∗ ȁ) for colloidal suspensions of different volume fractions (𝜙). 29 to check whether the TCS principle can be applied to a combination of both the 𝛼-relaxation and the 𝛽-relaxation processes in our model colloidal systems. For this purpose, 𝜙 = 0.40 was chosen as the reference volume fraction. The horizontal and the vertical shifts of the viscoelastic modulus values were performed to collapse them onto a single master curve. Specifically, the horizontal shift factors (𝑎𝜙 ) were first identified by constructing a master curve of the loss tangent (𝛿) as PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. and not to the 𝐺 ′ and the reduced loss modulus (𝐺 ′′ − 𝜔𝜂𝑠 ) values. In particular, it is of interest shown in Fig. 10. We note that a slope value of −0.5 in the low frequency region of the loss tangent master curve corresponds to the terminal region and the values at high volume fractions (𝜙 ≥ 0.45) were horizontally shifted with respect to this terminal slope line. Furthermore, for 𝜙 ≥ 0.45, an increase in the loss tangent values with an increase in the frequency indicates the emergence of the 𝛽-relaxation process. This exercise was then followed by the manual determination of the vertical shift factors (𝑏𝜙 ) that are required to collapse all of the 𝐺 ′ values at different volume fractions. The same shift factors were then used to construct a master curve of the 𝐺 ′′ values. Fig. 11 shows the master curves of the scaled elastic (𝑏𝜙 𝐺 ′ ) and the scaled loss (𝑏𝜙 𝐺 ′′ ) modulus as a function of the scaled frequency (𝑎𝜙 𝜔). As can be seen from the figure, the master curves obtained by the application of the TCS principle successfully extend the frequency range available to simulations, thus resulting in a linear viscoelastic (LVE) spectrum that spans nine decades of frequency in the case of 𝐺 ′ and more than four decades of frequency in the case of 𝐺 ′′ . The LVE spectrum clearly shows that the colloidal system response exhibits a transition from a fluid-like behavior (dominant 𝐺 ′′ ) in the terminal regime to more of a solid-like behavior (dominant 𝐺 ′ ) at the higher frequencies. However, the expected crossover region in 𝐺 ′′ (i.e. the region where a minimum in 𝐺 ′′ is observed) could not be obtained at the high volume fractions as it would require simulations at very low frequencies which are computationally expensive. The 30 region for the master curve, but the TCS principle is observed to work only up to a volume fraction of 𝜙 = 0.40. The loss (𝐺 ′′ ) modulus values fail to collapse onto the master curve at high volume fractions near the glass transition concentration (𝜙𝑔 ), due to the emergence and the overlap of the 𝛽-relaxation process with the 𝛼-relaxation process at these volume fractions, as was also observed by Peng et al. [76]. We do note that the 𝐺 ′ values are observed to seemingly superpose over the PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. expected terminal region scaling, i.e., 𝐺 ′ ~ 𝜔2 and 𝐺 ′′ ~ 𝜔1 , is observed in the low frequency entire frequency range. However, the 𝐺 ′ master curve obtained in this work should not be considered as a true master curve because the corresponding 𝐺 ′′ values do not superpose at high volume fractions. It is thus clear that only the 𝛼-relaxation process seems to follow the TCS principle, as was argued by Peng et al. [76]. The quality of the master curves obtained using the TCS principle can be further assessed by using the van Gurp-Palmen plot of the scaled complex modulus, 𝑏𝜙 ȁ𝐺 ∗ ȁ = 𝑏𝜙 √𝐺 ′ 2 + 𝐺 ′′ 2 [77] (see Fig 12). As can be seen from the figure, the phase lag values fall on a continuous curve from 𝜙 = 0.15 to 𝜙 = 0.40, thus indicating the validity of the TCS principle at low volume fractions. Furthermore, the phase lag values for 𝜙 > 0.40 do not follow a continuous trend, a finding that is consistent with the observation that the TCS principle fails at high volume fractions near 𝜙𝑔 . A phase lag value of 𝛿 ≈ 1.57 𝑟𝑎𝑑 ≈ 90° corresponding to the terminal relaxation is also obtained, which is consistent with the terminal slopes seen in the master curves of the LVE spectrum, where the TCS principle is valid. To further validate the master curves, the Kramers-Kronig relations [78], as given by Eqs. (8) and (9), were used to predict the 𝐺 ′ values from the 𝐺 ′′ values and vice-versa. 𝐺 ′ (𝜔) = 𝐺 ′ (0) − 2𝜔2 𝜋 ∞ 𝐺 ′′ (𝑢)⁄𝑢−𝐺 ′′ (𝜔)⁄𝜔 𝒫 ∫0 31 𝑢2 −𝜔2 𝑑𝑢 (8) This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 (𝒂) (𝒃) FIG. 𝟏𝟑. Comparison of the Kramers-Kronig predictions (yellow plus-squares) with the master curves of (𝑎) the scaled elastic (𝑏𝜙 𝐺 ′ ) and (𝑏) the scaled loss (𝑏𝜙 𝐺 ′′ ) modulus values. 32 2𝜔 𝜋 ∞ 𝐺 ′ (𝑢)−𝐺 ′ (𝜔) 𝒫 ∫0 𝑢2 −𝜔2 𝑑𝑢 (9) where 𝒫 is the Cauchy principal value. Fig. 13 compares the Kramers-Kronig predictions of the scaled elastic (𝑏𝜙 𝐺 ′ ) and the scaled loss (𝑏𝜙 𝐺 ′′ ) modulus values with the master curves obtained using the TCS principle. We note that the Kramers-Kronig predictions of 𝑏𝜙 𝐺 ′′ were obtained only in the low frequency region where the TCS principle is successfully applicable. As can be PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 𝐺 ′′ (𝜔) = seen in Figs. 13(𝑎) and 13(𝑏), the 𝑏𝜙 𝐺 ′ and the 𝑏𝜙 𝐺 ′′ predictions obtained using Eqs. (8) and (9) show small variation from the values obtained from the master curves of 𝑏𝜙 𝐺 ′ and 𝑏𝜙 𝐺 ′′ at all frequencies. We attribute the differences in the Kramers-Kronig predictions of 𝑏𝜙 𝐺 ′ to the absence of 𝑏𝜙 𝐺 ′′ values at high frequencies whereas the differences in 𝑏𝜙 𝐺 ′′ values are attributed to the very large uncertainties in 𝑏𝜙 𝐺 ′ values in the low frequency region, which make it difficult to obtain the 𝑏𝜙 𝐺 ′′ values accurately. In summary, the TCS principle is shown to be applicable to the simulation results for the viscoelastic moduli of the model colloidal systems only for the volume fractions well below the glass transition concentration (i.e., for the terminal region related to the 𝛼-relaxation process). 1. Physical significance of the shift factors The shift factors can be related to the relative zero shear viscosity of the colloidal suspensions as follows [79]: 𝑎𝜙 𝑏𝜙 = 𝜂𝑟,𝜙 𝜂𝑟,𝑟𝑒𝑓 (10) where 𝜂𝑟,𝜙 and 𝜂𝑟,𝑟𝑒𝑓 are the relative zero shear viscosity values at the volume fraction of interest and the reference volume fraction, respectively. Fig. 14 shows a comparison of the plots of the 33 PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 𝑎 FIG. 𝟏𝟒. Comparison of the ratio of shift factors ( 𝜙 ) and relative zero shear viscosity 𝑏 𝜂𝑟,𝜙 (𝜂 𝑟,𝑟𝑒𝑓 𝜙 ) as a function of volume fraction (𝜙). The dashed violet line is the VFT fit to the ratio of shift factors. Inset: Horizontal (𝑎𝜙 ) and vertical (𝑏𝜙 ) shift factors as a function of volume fraction (𝜙). Also, the horizontal shift factors (𝑎𝜙,𝛽 ) that were used to superpose the 𝐺 ′′ values in the 𝛽-relaxation region are shown as pink cross-circles. The dashed violet line is the VFT fit to the 𝑎𝜙 values. 34 𝑏𝜙 𝑟,𝑟𝑒𝑓 ) as a function of the colloidal volume fraction (𝜙). It is seen that, within the statistical uncertainties, values of the ratio of shift factors are in good agreement with the relative zero shear viscosity values up to 𝜙 = 0.40 but they start to deviate from the relative zero shear viscosity values at 𝜙 = 0.45. We note that the TCS principle was also found to be applicable only up to the volume fraction of 𝜙 = 0.40. Also, the PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 𝜂𝑟,𝜙 𝑎𝜙 ratio of shift factors ( ) and the relative zero shear viscosity ( 𝜂 horizontal shift factors are observed to increase with an increase in the volume fraction up to 𝜙 = 0.49 whereas the vertical shift factors increase up to = 0.40 and attain almost a constant value at higher volume fractions (see inset of Fig. 14). Furthermore, the TCS principle was successfully applied to superpose the 𝛽-relaxation process individually and the horizontal shift factors (𝑎𝜙,𝛽 ) corresponding to this superposition are also shown in the inset of Fig. 14. As can be seen, these 𝑎𝜙,𝛽 values do not show a dependence on the volume fraction as was also seen in the experimental and simulation study by Peng et al. [76]. The volume fraction dependence of the ratio of shift factors or the relaxation time can be described [76, 80] by the Vogel-Fulcher-Tammann (VFT) equation as follows [81-83]: 𝑎𝜙 𝑏𝜙 = 𝑒𝑥𝑝 [𝐵 + 𝐶 𝜙∞ −𝜙 ] (11) where B and C are the VFT constants, and 𝜙∞ represents the volume fraction at which the motion of the colloidal particles freezes completely (i.e., the random close packing, RCP). We note that a good fit of the William-Landel-Ferry (WLF) equation [84, 85] (which is mathematically equivalent to the VFT equation) to the shift factors of the unjammed systems was reported by Wen et al. [75] in the application of TCS principle to study the dynamics of soft colloidal glasses. As can be seen from the Fig. 14, the VFT equation provides an excellent fit to the ratio of the shift 35 fractions of 𝜙 < 0.455). The best fit parameter values of 𝐵 = −3.886 and 𝐶 = 0.483 were obtained from this VFT fit to the 𝑎𝜙 𝑏𝜙 ratio. These values are similar to those obtained in a recent study of colloidal glasses by one of us [76]. Furthermore, a fit value of 𝜙∞ = 0.524 was obtained which is slightly higher than the glass transition concentration (𝜙𝑔 ) value of 0.519 that is obtained PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. factors at volume fractions lower than the glass transition concentration (i.e., unjammed volume using the MCT fit to the 𝜂𝑟 values (see Sec. III B). We note that the VFT equation can also be used to represent the horizontal shift factors (𝑎𝜙 ) used in the construction of the master curves. Such an attempt to fit the VFT equation to the 𝑎𝜙 values resulted in the parameter values of 𝐵 = −12.17, 𝐶 = 3.357, and 𝜙∞ = 0.676 (see inset of Fig. 14). The 𝜙∞ value obtained from this fit is much higher than expected and is not consistent with the previous observations made for our model systems. Also, the VFT equation could not properly represent the 𝑎𝜙 values at low volume fractions. These findings signify the importance of accounting for the vertical shift factors (𝑏𝜙 ) (i.e., using the 𝑎𝜙 𝑏𝜙 ratio instead of just 𝑎𝜙 ) for the purpose of accurately quantifying the viscosity behavior as a function of volume fraction. Both the 𝑎𝜙 and the 𝑏𝜙 values contain information related to the relaxation modes of the model colloidal systems and both must be taken into account whenever the vertical shift factors along with the horizontal shift factors are used in the construction of the master curves of the elastic and the loss modulus values using the TCS principle. IV. SUMMARY AND CONCLUSIONS The interparticle interactions play a significant role in governing the rheology of nanocolloidal suspensions due to the comparable sizes of the solvent and the colloidal particles in these systems. 36 solvent and solvent-solvent interactions. The model is used to investigate the rheology of nanocolloidal suspensions using MD simulations. The rheological features of our nanocolloidal model systems with a 5: 1 colloid-solvent size ratio were found to be consistent with those observed in the systems with a large size difference between the solvent and the colloidal particles. Specifically, for 𝜙 ≤ 0.45, the relative zero shear viscosity (𝜂𝑟 ) values of these systems were observed to be in good agreement with those reported in the previous macroscopic experimental PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. In this work, we have used a particulate molecular model that explicitly accounts for these colloid- and mesoscopic simulation studies. The elastic (𝐺 ′ ) and the loss (𝐺 ′′ ) modulus values of these model colloidal suspensions showed the expected transition from a fluid-like behavior at low volume fractions to more of a solid-like behavior at high volume fractions. Molecular simulations can treat the interparticle interactions via the use of explicit solvent models but because of this aspect, the simulations are restricted to small time or frequency scale windows. Our results showed that similar to experimental rheology, the frequency scale accessible to simulations can be extended by applying the time-concentration superposition (TCS) principle to the colloid rheological properties but the TCS principle works only for the 𝛼-relaxation process. Application of the TCS principle resulted in a linear viscoelastic (LVE) spectrum covering nine decades of frequency for the elastic modulus and more than four decades of frequency for the loss modulus values. However, the TCS principle failed at high volume fractions (𝜙 > 0.40) near the glass transition concentration (𝜙𝑔 ) regime, as evident from the master curve of the loss modulus. The validity of the TCS principle was judged by using several tests: The van Gurp-Palmen plot showed a continuous behavior of the phase lag as a function of the complex modulus and the master curves followed the Kramers-Kronig predictions. Furthermore, the volume fraction dependence of the ratio of shift factors followed the VFT equation and was also in quantitative 37 particulate molecular model of the nanocolloidal suspensions captures the known rheological behavior of colloidal systems in terms of viscosity and viscoelasticity and the time-concentration superposition (TCS) principle can be used as a viable tool to overcome the challenges posed by the accessibility to a limited timescale in simulations. The explicit solvent molecular model provides the ability to study the fundamental physics PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. agreement with the relative zero shear viscosity values. In summary, our results show that a fully related questions in probe rheology that result when the probe particle size and the suspension medium particle size become comparable. Such an effort that combines probe rheology simulations [86, 87] with the molecular models of nanocolloidal suspensions is currently underway in our group. The interactions between different types of particles were represented using the basic, purely repulsive WCA potential in this work. The ability of the model to handle the effect of specific interactions was demonstrated by making a simple change in the interaction potential from the WCA to the LJ potential (which also has attractive interactions); this change resulted in approximately a factor of two increase in the viscosity of the model system. As shown in the previous literature studies [24, 48], the local structure of the colloidal systems (e.g. aggregation/dispersion) can be tuned by changing the length scale of the different interactions in the systems. Furthermore, the model used in this work can be readily extended by incorporating the appropriate intermolecular interactions (e.g. electrostatic interactions) between the particles in the case if the interactions in the system are strong and thus one is interested in studying the effects of specific chemistry on the structure as well as the rheology of a particular nanocolloidal suspension. Such a study can be carried out using atomistically detailed models of the nanocolloidal systems. However, this ability to handle specific interactions and short range structure comes with a potentially severe increase in the computational cost of the simulations 38 representation of the solvent. The explicit solvent model that is employed here will only be useful when specific chemical interactions and nanoscale molecular structure play an important role in governing the system rheology thus rendering continuum treatment of the solvent inadequate. For complex fluid systems where both nanoscale and mesoscale effects are important, the results from these explicit solvent model simulations could serve as an input to the mesoscale simulation techniques such as BD or SD simulations. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. compared to the cost of mesoscopic simulation techniques such as BD and SD that use a continuum 39 This material is based upon work partially supported by the National Science Foundation under the grant NSF DMR-1611328. The computational resources provided by the High Performance Computing Center (HPCC) at Texas Tech University and Texas Advanced Computing Center (TACC) at The University of Texas at Austin are acknowledged. This work also used the Extreme PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. ACKNOWLEDGMENTS Science and Engineering Discovery Environment (XSEDE) allocation time on Stampede-2 at TACC, which is supported by the National Science Foundation under the grant NSF ACI1548562. 40 1. Zia, R.N., "Active and passive microrheology: Theory and simulation," Annual Review of Fluid Mechanics 50, 371-405 (2018). 2. Pusey, P.N., and W. van Megen, "Phase behaviour of concentrated suspensions of nearly hard colloidal spheres," Nature 320, 340 (1986). 3. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. REFERENCES Zia, R.N., and J.F. Brady, "Microviscosity, microdiffusivity, and normal stresses in colloidal dispersions," Journal of Rheology 56, 1175-1208 (2012). 4. Swan, J.W., R.N. Zia, and J.F. Brady, "Large amplitude oscillatory microrheology," Journal of Rheology 58, 1-41 (2014). 5. Khair, A.S., and J.F. Brady, "“Microviscoelasticity” of colloidal dispersions," Journal of Rheology 49, 1449-1481 (2005). 6. Squires, T.M., and J.F. Brady, "A simple paradigm for active and nonlinear microrheology," Physics of Fluids 17, 073101 (2005). 7. Cichocki, B., and B.U. Felderhof, "Linear viscoelasticity of colloidal suspensions," Physical Review A 46, 7723-7732 (1992). 8. Varga, Z., and J.W. Swan, "Linear viscoelasticity of attractive colloidal dispersions," Journal of Rheology 59, 1271-1298 (2015). 9. Brady, J.F., A.S. Khair, and M. Swaroop, "On the bulk viscosity of suspensions," Journal of Fluid Mechanics 554, 109-123 (2006). 10. Nazockdast, E., and J.F. Morris, "Microstructural theory and the rheology of concentrated colloidal suspensions," Journal of Fluid Mechanics 713, 420-452 (2012). 11. Mason, T.G., and D.A. Weitz, "Linear viscoelasticity of colloidal hard sphere suspensions near the glass transition," Physical Review Letters 75, 2770-2773 (1995). 41 Sriram, I., A. Meyer, and E.M. Furst, "Active microrheology of a colloidal suspension in the direct collision limit," Physics of Fluids 22, 062003 (2010). 13. Habdas, P., D. Schaar, A.C. Levitt, and E.R. Weeks, "Forced motion of a probe particle near the colloidal glass transition," EPL (Europhysics Letters) 67, 477 (2004). 14. Lee, M., M. Alcoutlabi, J.J. Magda, C. Dibble, M.J. Solomon, X. Shi, and G.B. McKenna, "The effect of the shear-thickening transition of model colloidal spheres on the sign of N1 and on the radial pressure profile in torsional shear flows," Journal of Rheology 50, 293- PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 12. 311 (2006). 15. Bender, J., and N.J. Wagner, "Reversible shear thickening in monodisperse and bidisperse colloidal dispersions," Journal of Rheology 40, 899-916 (1996). 16. Lee, Y.S., and N.J. Wagner, "Dynamic properties of shear thickening colloidal suspensions," Rheologica Acta 42, 199-208 (2003). 17. Meeker, S.P., W.C.K. Poon, and P.N. Pusey, "Concentration dependence of the low-shear viscosity of suspensions of hard-sphere colloids," Physical Review E 55, 5718-5722 (1997). 18. Wilson, L.G., A.W. Harrison, A.B. Schofield, J. Arlt, and W.C.K. Poon, "Passive and active microrheology of hard-sphere colloids," The Journal of Physical Chemistry B 113, 3806-3812 (2009). 19. Seth, J.R., L. Mohan, C. Locatelli-Champagne, M. Cloitre, and R.T. Bonnecaze, "A micromechanical model to predict the flow of soft particle glasses," Nature Materials 10, 838 (2011). 20. Ermak, D.L., and J.A. McCammon, "Brownian dynamics with hydrodynamic interactions," The Journal of Chemical Physics 69, 1352-1360 (1978). 42 Carpen, I.C., and J.F. Brady, "Microrheology of colloidal dispersions by Brownian dynamics simulations," Journal of Rheology 49, 1483-1502 (2005). 22. Brady, J.F., and G. Bossis, "Stokesian dynamics," Annual Review of Fluid Mechanics 20, 111-157 (1988). 23. Nazockdast, E., and J.F. Morris, "Effect of repulsive interactions on structure and rheology of sheared colloidal dispersions," Soft Matter 8, 4223-4234 (2012). 24. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 21. Padding, J.T., and A.A. Louis, "Hydrodynamic interactions and Brownian forces in colloidal suspensions: Coarse-graining over time and length scales," Physical Review E 74, 031402 (2006). 25. Pryamitsyn, V., and V. Ganesan, "A coarse-grained explicit solvent simulation of rheology of colloidal suspensions," The Journal of Chemical Physics 122, 104906 (2005). 26. Malevanets, A., and R. Kapral, "Solute molecular dynamics in a mesoscale solvent," The Journal of Chemical Physics 112, 7260-7269 (2000). 27. Khare, R., and D. Sundaravadivelu Devarajan, "Molecular simulations of nanocolloids," Current Opinion in Chemical Engineering 16, 86-91 (2017). 28. Gossard, A., F. Frances, and C. Aloin, "Rheological properties of TiO2 suspensions varied by shifting the electrostatic inter-particle interactions with an organic co-solvent," Colloids and Surfaces A: Physicochemical and Engineering Aspects 522, 425-432 (2017). 29. Wang, T., M. Ni, Z. Luo, C. Shou, and K. Cen, "Viscosity and aggregation structure of nanocolloidal dispersions," Chinese Science Bulletin 57, 3644-3651 (2012). 30. Pamies, R., C. Espejo, F.J. Carrión, A. Morina, A. Neville, and M.D. Bermúdez, "Rheological behavior of multiwalled carbon nanotube-imidazolium tosylate ionic liquid dispersions," Journal of Rheology 61, 279-289 (2017). 43 Tan, Y., Y. Song, and Q. Zheng, "Hydrogen bonding-driven rheological modulation of chemically reduced graphene oxide/poly(vinyl alcohol) suspensions and its application in electrospinning," Nanoscale 4, 6997-7005 (2012). 32. Shim, Y.H., K.E. Lee, T.J. Shin, S.O. Kim, and S.Y. Kim, "Tailored colloidal stability and rheological properties of graphene oxide liquid crystals with polymer-induced depletion attractions," ACS Nano 12, 11399-11406 (2018). 33. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 31. Chau, M., S.E. Sriskandha, D. Pichugin, H. Thérien-Aubin, D. Nykypanchuk, G. Chauve, M. Méthot, J. Bouchard, O. Gang, and E. Kumacheva, "Ion-mediated gelation of aqueous suspensions of cellulose nanocrystals," Biomacromolecules 16, 2455-2462 (2015). 34. Saluja, A., A.V. Badkar, D.L. Zeng, S. Nema, and D.S. Kalonia, "Ultrasonic storage modulus as a novel parameter for analyzing protein-protein interactions in high protein concentration solutions: Correlation with static and dynamic light scattering measurements," Biophysical Journal 92, 234-244 (2007). 35. Sun, S.m., Y. Song, and Q. Zheng, "pH-induced rheological changes for semi-dilute solutions of wheat gliadins," Food Hydrocolloids 22, 1090-1096 (2008). 36. Woldeyes, M.A., L.L. Josephson, D.L. Leiske, W.J. Galush, C.J. Roberts, and E.M. Furst, "Viscosities and protein interactions of bispecific antibodies and their monospecific mixtures," Molecular Pharmaceutics 15, 4745-4755 (2018). 37. Sharma, A.K., A.K. Tiwari, and A.R. Dixit, "Rheological behaviour of nanofluids: A review," Renewable and Sustainable Energy Reviews 53, 779-791 (2016). 38. Lin, S., and D. Blankschtein, "Role of the bile salt surfactant sodium cholate in enhancing the aqueous dispersion stability of single-walled carbon nanotubes: A molecular dynamics simulation study," The Journal of Physical Chemistry B 114, 15616-15625 (2010). 44 Bari, R., D. Parviz, F. Khabaz, C.D. Klaassen, S.D. Metzler, M.J. Hansen, R. Khare, and M.J. Green, "Liquid phase exfoliation and crumpling of inorganic nanosheets," Physical Chemistry Chemical Physics 17, 9383-9393 (2015). 40. Habib, T., D. Sundaravadivelu Devarajan, F. Khabaz, D. Parviz, T.C. Achee, R. Khare, and M.J. Green, "Cosolvents as liquid surfactants for boron nitride nanosheet (BNNS) dispersions," Langmuir 32, 11591-11599 (2016). 41. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 39. Heine, D.R., M.K. Petersen, and G.S. Grest, "Effect of particle shape and charge on bulk rheology of nanoparticle suspensions," The Journal of Chemical Physics 132, 184509 (2010). 42. in ’t Veld, P.J., M.K. Petersen, and G.S. Grest, "Shear thinning of nanoparticle suspensions," Physical Review E 79, 021401 (2009). 43. Weeks, J.D., D. Chandler, and H.C. Andersen, "Role of repulsive forces in determining the equilibrium structure of simple liquids," The Journal of Chemical Physics 54, 5237-5247 (1971). 44. Abascal, J.L.F., and C. Vega, "A general purpose model for the condensed phases of water: TIP4P/2005," The Journal of Chemical Physics 123, 234505 (2005). 45. Periole, X., and S.-J. Marrink, The Martini Coarse-Grained Force Field, in Biomolecular Simulations: Methods and Protocols, L. Monticelli and E. Salonen, Editors. 2013, Humana Press: Totowa, NJ. p. 533-565. 46. Schmidt, J.R., and J.L. Skinner, "Hydrodynamic boundary conditions, the Stokes–Einstein law, and long-time tails in the Brownian limit," The Journal of Chemical Physics 119, 8062-8068 (2003). 45 Kohale, S.C., and R. Khare, "Molecular simulation of cooperative hydrodynamic effects in motion of a periodic array of spheres between parallel walls," The Journal of Chemical Physics 129, 164706 (2008). 48. Tomilov, A., A. Videcoq, T. Chartier, T. Ala-Nissilä, and I. Vattulainen, "Tracer diffusion in colloidal suspensions under dilute and crowded conditions with hydrodynamic interactions," The Journal of Chemical Physics 137, 014503 (2012). 49. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 47. Khabaz, F., and R. Khare, "Molecular simulations of asphalt rheology: Application of time–temperature superposition principle," Journal of Rheology 62, 941-954 (2018). 50. Evans, D.J., and G. Morriss, Statistical Mechanics of Nonequilibrium Liquids (Cambridge University Press, 2008). 51. Cornelis Elisa Struik, L., Physical Aging In Amorphous Polymers and Other Materials 1978). 52. Peng, X., and G.B. McKenna, "Physical aging and structural recovery in a colloidal glass subjected to volume-fraction jump conditions," Physical Review E 93, 042603 (2016). 53. Shinoda, W., M. Shiga, and M. Mikami, "Rapid estimation of elastic constants by molecular dynamics simulation under constant stress," Physical Review B 69, 134103 (2004). 54. Plimpton, S., "Fast parallel algorithms for short-range molecular dynamics," Journal of Computational Physics 117, 1-19 (1995). 55. Larson, R.G., The Structure and Rheology of Complex Fluids (Oxford University, 1998). 56. Grest, G.S., Q. Wang, P.i.t. Veld, and D.J. Keffer, "Effective potentials between nanoparticles in suspension," The Journal of Chemical Physics 134, 144902 (2011). 46 Khabaz, F., T. Liu, M. Cloitre, and R.T. Bonnecaze, "Shear-induced ordering and crystallization of jammed suspensions of soft particles glasses," Physical Review Fluids 2, 093301 (2017). 58. Kohale, S.C., and R. Khare, "Molecular dynamics simulation study of friction force and torque on a rough spherical particle," The Journal of Chemical Physics 132, 234706 (2010). 59. Russel, W.B., N.J. Wagner, and J. Mewis, "Divergence in the low shear viscosity for Brownian hard-sphere dispersions: At random close packing or the glass transition?," PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 57. Journal of Rheology 57, 1555-1567 (2013). 60. Sriram, I., E.M. Furst, R.J. DePuit, and T.M. Squires, "Small amplitude active oscillatory microrheology of a colloidal suspension," Journal of Rheology 53, 357-381 (2009). 61. Segrè, P.N., S.P. Meeker, P.N. Pusey, and W.C.K. Poon, "Viscosity and structural relaxation in suspensions of hard-sphere colloids," Physical Review Letters 75, 958-961 (1995). 62. Cheng, Z., J. Zhu, P.M. Chaikin, S.-E. Phan, and W.B. Russel, "Nature of the divergence in low shear viscosity of colloidal hard-sphere dispersions," Physical Review E 65, 041405 (2002). 63. Shikata, T., and D.S. Pearson, "Viscoelastic behavior of concentrated spherical suspensions," Journal of Rheology 38, 601-616 (1994). 64. Foss, D.R., and J.F. Brady, "Structure, diffusion and rheology of Brownian suspensions by Stokesian dynamics simulation," Journal of Fluid Mechanics 407, 167-200 (2000). 65. Phung, T.N., "Behavior of concentrated colloidal suspensions by Stokesian dynamics simulation," PhD Thesis, California Institute of Technology, (1993). 47 Pan, W., B. Caswell, and G.E. Karniadakis, "Rheology, microstructure and migration in Brownian colloidal suspensions," Langmuir 26, 133-142 (2010). 67. Krieger, I.M., and T.J. Dougherty, "A mechanism for non-Newtonian flow in suspensions of rigid spheres," Transactions of the Society of Rheology 3, 137-152 (1959). 68. Lionberger, R.A., and W.B. Russel, "Microscopic theories of the rheology of stable colloidal dispersions," Advances in Chemical Physics 111, 399 (2000). 69. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 66. Bengtzelius, U., W. Gotze, and A. Sjolander, "Dynamics of supercooled liquids and the glass transition," Journal of Physics C: Solid State Physics 17, 5915-5934 (1984). 70. Al-Hadithi, T.S.R., H.A. Barnes, and K. Walters, "The relationship between the linear (oscillatory) and nonlinear (steady-state) flow properties of a series of polymer and colloidal systems," Colloid and Polymer Science 270, 40-46 (1992). 71. Bertsch, P., A. Sánchez-Ferrer, M. Bagnani, S. Isabettini, J. Kohlbrecher, R. Mezzenga, and P. Fischer, "Ion-induced formation of nanocrystalline cellulose colloidal glasses containing nematic domains," Langmuir 35, 4117-4124 (2019). 72. Cox, W.P., and E.H. Merz, "Correlation of dynamic and steady flow viscosities," Journal of Polymer Science 28, 619-622 (1958). 73. Brady, J.F., "The rheological behavior of concentrated colloidal dispersions," The Journal of Chemical Physics 99, 567-581 (1993). 74. Cicuta, P., E.J. Stancik, and G.G. Fuller, "Shearing or compressing a soft glass in 2D: Time-concentration superposition," Physical Review Letters 90, 236101 (2003). 75. Wen, Y.H., J.L. Schaefer, and L.A. Archer, "Dynamics and rheology of soft colloidal glasses," ACS Macro Letters 4, 119-123 (2015). 48 Peng, X., J.G. Wang, Q. Li, D. Chen, R.N. Zia, and G.B. McKenna, "Exploring the validity of time-concentration superposition in glassy colloids: Experiments and simulations," Physical Review E 98, 062602 (2018). 77. Qian, Z., and G.B. McKenna, "Expanding the application of the van Gurp-Palmen plot: New insights into polymer melt rheology," Polymer 155, 208-217 (2018). 78. Booij, H.C., and G.P.J.M. Thoone, "Generalization of Kramers-Kronig transforms and some approximations of relations between viscoelastic quantities," Rheologica Acta 21, PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 76. 15-24 (1982). 79. Markovitz, H., "Superposition in rheology," Journal of Polymer Science: Polymer Symposia 50, 431-456 (1975). 80. Li, Q., X. Peng, and G.B. McKenna, "Long-term aging behaviors in a model soft colloidal system," Soft Matter 13, 1396-1404 (2017). 81. Vogel, H., "Das Temperaaturabhä ngigkeits gesetz der Viskositä t Flüssigkeiten," Phys. Zeit. 22, 645-646 (1921). 82. Fulcher, G.S., "Analysis of recent measurements of the viscosity of glasses," Journal of the American Ceramic Society 8, 339-355 (1925). 83. Tammann, G., "Glasses as supercooled liquids," Journal of the Society of Glass Technology. 9, 166-185 (1925). 84. Williams, M.L., R.F. Landel, and J.D. Ferry, "The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids," Journal of the American Chemical Society 77, 3701-3707 (1955). 85. Rubinstein, M., and R.H. Colby, Polymer Physics (Oxford University, 2003). 49 Karim, M., T. Indei, J.D. Schieber, and R. Khare, "Determination of linear viscoelastic properties of an entangled polymer melt by probe rheology simulations," Physical Review E 93, 012501 (2016). 87. Karim, M., S.C. Kohale, T. Indei, J.D. Schieber, and R. Khare, "Determination of viscoelastic properties by analysis of probe-particle motion in molecular simulations," Physical Review E 86, 051501 (2012). PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 86. 50 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142 This is the author’s peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1122/1.5125142