Multiscale Molecular Simulations of Polymer-Matrix Nanocomposites

Following the substantial progress in molecular simulations of polymer-matrix nanocomposites, now is the time to reconsider this topic from a critical point of view. A comprehensive survey is reported herein providing an overview of classical molecular simulations, reviewing their major achievements in modeling polymer matrix nanocomposites, and identifying several open challenges. Molecular simulations at multiple length and time scales, working hand-in-hand with sensitive experiments, have enhanced our understanding of how nanofillers alter the structure, dynamics, thermodynamics, rheology and mechanical properties of the surrounding polymer matrices.


Introduction
Polymer-matrix nanocomposites (PNCs) have drawn intense research interest over the last decade owing to both the rich fundamental physics associated with mixing macromolecules and particles and their unique mechanical, optical, magnetic and other material properties [1]. Driven by the need to develop functionally superior materials, significant effort has been invested in understanding than 10 −3 within a polymer matrix. They are thus characterized by particle number densities n = 3 ∕ 4 3 p ≈ 10 20 m −3 , interfacial areas per unit volume 3 ∕ p ≈ 10 6 m −1 , and interparticle spacings, −1∕3 n − 2 p ≈ 100 nm that are commensurate with the particle dimensions, p and the radii of gyration of matrix chains, R g ≈ 10 nm.
The practice of adding nanoscale filler particles to reinforce polymeric materials can be traced back to the early years of the composite industry, in the second half of the 19th century. Charles Goodyear, inventor of vulcanized rubber, attempted to prepare nanoparticle-toughened automobile tires by blending carbon black, zinc oxide, and/or magnesium sulfate particles with vulcanized rubber [20]. Another example was the clay-reinforced resin known as Bakelite that was introduced in the early 1900s as one of the first mass-produced polymer-nanoparticle composites and fundamentally transformed the nature of practical household materials [21][22][23][24]. Then, a long period of time passed till the early 1990s when it was first demonstrated that the thermal and mechanical properties of Nylon-6 were improved by the addition of a few percent (2-4 % w/w) mica-type layered silicates to the extent that it could be used in an automotive engine compartment [25,26].
Even though some property improvements have been achieved in nanocomposites, nanoparticle dispersion is difficult to control, with both thermodynamic and kinetic processes playing significant roles. It has been demonstrated that dispersed spherical nanoparticles can yield a range of multifunctional behavior, including a viscosity decrease, reduction of thermal degradation, increased mechanical damping, enriched electrical and/or magnetic performance, and control of thermomechanical properties [27][28][29][30][31]. The tailor-made properties of these systems are very important to the manufacturing procedure, as they fully overcome many of the existing operational limitations. As a final product, a polymeric matrix enriched with dispersed particles may have better properties than the neat polymeric material and can be used in more demanding and novel applications. Therefore, an understanding and quantitative description of the physicochemical properties of these materials is of major importance for their successful production.
As part of this renewed interest in nanocomposites, researchers also began seeking design rules that would allow them to engineer materials that combine the desirable properties of nanoparticles and polymers. In light of the diversity of polymers and nanoparticles, the potential for use of PNCs is nearly limitless. The ensuing research revealed a number of key challenges in producing nanocomposites that exhibit a desired behavior. The greatest stumbling block to the large-scale production and commercialization of nanocomposites is the dearth of cost-effective methods for controlling the dispersion of the nanoparticles in polymeric hosts. The nanoscale particles typically aggregate, which negates any benefits associated with the nanoscopic dimension. PNCs generally possess nonequilibrium morphologies due to the complex interplay of enthalpic and entropic interactions leading to particle aggregation, particle bridging interactions, and phase separation at various length scales [32,33]. The second challenge is associated with understanding and predicting property enhancements in these materials, which are intimately connected to their morphology.
Nanocomposite research has recently expanded to consider more complicated systems involving polymer blends and block copolymers, where novel electrical, magnetic and optical properties arise [15,34,35].

Multiscale Modeling
Understanding the fascinating and complex structure and dynamics of polymeric materials has been an ongoing challenge for many decades. From the point of view of molecular simulations, the spectrum of length and time scales associated with polymer melts of long chains poses a formidable challenge to studying their long-time dynamics [36,37]. The topological constraints arising from chain connectivity and uncrossability (entanglements) dominate intermediate and long-time relaxation [38] and transport phenomena when polymers become sufficiently long. Atomistic molecular simulations of dense phases of soft matter prove to be difficult for many systems across length and time scales of practical interest. Even coarse-grained particle-based simulation methods may not be applicable due to the lack of faithful descriptions of polymer-polymer and polymer-surface interactions. Since complex interactions between constituent phases at the atomic level ultimately manifest themselves in macroscopic properties, a broad range of length and time scales must be addressed and a combination of modeling techniques is therefore required to simulate meaningfully the bulk-level behavior of nanocomposites [9].
Soft condensed matter is a relatively new term describing a huge class of rather different materials such as colloids, polymers, membranes, complex molecular assemblies, complex fluids etc. Though these materials are rather different in their structures, there is one unifying aspect, which makes it very reasonable to treat such systems from a common point of view. Compared to "hard matter" the characteristic energy density is much smaller. While the typical energy of a chemical bond (C-C bond) is about 10 −18 J ≈ 250k B T at room temperature of 300 K, the nonbonded interactions are of the order of k B T and allow for strong density fluctuations even though the molecular connectivity is never affected (k B is the Boltzmann's constant). It is instructive to compare the cohesive energy density, which gives a first estimate of the elastic constants, between a typical "hard matter" crystal to soft matter. The ratio between the two shows that polymeric systems are typically 100-10,000 times softer than classical crystals. As a consequence the average thermal energyk B T is not negligible for these systems any more, but rather defines the essential energy scale. This means that entropy, which typically contributes to the free energy a term of the order of k B T per degree of freedom, plays a crucial role. Especially in the case of macromolecules, this is mainly intramolecular entropy, which for a linear polymer of length N contributes to the free energy a term of order Nk B T, representing about 90% of the free energy of polymeric materials [39]. As an immediate consequence it is clear that typical quantum chemical electronic structure calculations (Hartree-Fock or DFT) which focus on obtaining the energy as a function of nuclear coordinates cannot be sufficient to characterize soft condensed matter and will even be less sufficient to properly predict/interpret macroscopic properties. Molecular theoretical and simulation methods which incorporate entropic effects are required for this.
The length and time scales governing polymer physics range from Å and femtoseconds for the vibrations of atomic bonds to millimeters and seconds for crack propagation in polymer composites. The entities used as basic degrees of freedom are: electrons (quantum chemistry), atoms (classical forcefields), monomers or groups of monomers (coarsegrained or mesoscopic models) and entire polymer chains (soft fluids). All these methods and many others have been applied side by side to polymers. Until recently, however, multiscale methods with rigorous bridging between the different scales have been few.

Atomistic Molecular Dynamics (MD)
The stepping stone of classical molecular simulations is atomistic Molecular Dynamics (MD). As accurate MD potentials are developed for a broad range of materials based on quantum chemistry calculations and with the increase of supercomputer performance, atomistic MD simulations have become a very powerful tool for analyzing complex physical phenomena in polymeric materials, including dynamics, viscosity and shear thinning. However, as discussed above, entangled polymer systems are characterized by a wide range of spatial and temporal scales. It is still not feasible to equilibrate atomistic MD simulations of highly entangled polymer chain systems, due to their long relaxation times, long-range electrostatic interactions and tremendous number of atoms. The atomistic MD model for such a system, with a typical size of about a micrometer and a relaxation time on the scale of microseconds (or even up to the scale of seconds for long-chain polymer melts), would consist of billions of atoms and would require billions of time steps to run, which is obviously beyond the capability of the technique, even with the most sophisticated supercomputers available today.

Monte Carlo (MC)
A robust sampling of the configuration space of polymeric substances is a prerequisite for the reliable prediction of their physical properties. The constraints posed by atomistic MD simulations can be overcome by resorting to MC simulations, which enable us to use the complete arsenal of equilibrium statistical mechanics, e.g. perform sampling in all sorts of ensembles [36,37,[40][41][42]. Through the design of efficient unphysical moves, configurational sampling can be dramatically enhanced. MC moves such as concerted rotation [43], configurational bias [44,45], and internal configurational bias [46] have thus successfully addressed the problem of equilibrating polymer systems of moderate chain lengths.
Even these moves prove incapable of providing equilibration when applied to long-chain polymer melts, however. A solution to this problem was given by the development and efficient implementation of a chain connectivity-altering MC move, end-bridging [47,48]. Using end-bridging, atomistic systems consisting of a large number of long chains, up to C 6000 , have been simulated in full atomistic detail [48,49]. Despite its efficiency in equilibrating long-chain polymer melts, end-bridging cannot equilibrate monodisperse polymer melts; a finite degree of polydispersity is necessary for the move to operate. While this is not a drawback in modeling industrial polymers, which are typically polydisperse, an ability to equilibrate strictly monodisperse polymers is highly desirable for comparing against theory or model experimental systems. Morover, end-bridginig relies on the existence of chain ends, rendering itself inappropriate for dense phases of chains with nonlinear architectures. These limitations have been overcome by the introduction of Double Bridging (DB) and Intramolecular Double Rebridging (IDR) [50,51]. The key innovation of those moves is the construction of two bridging trimers between two different chains, as far as the former is concerned, or along the same chain, as far as the latter move is concerned, thus preserving the initial chain lengths.
MC simulations using atomistic forcefields have inherent limitations, as Doxastakis et al. have shown [52]. The hard interactions between atoms reduce the acceptance rate of the moves. Thus, it is essential to resort to parallel tempering techniques in order to allow motion of the system in its phase space [53].

Coarse Graining (CG)
Polymers show a hierarchy of length and time scales. However, the connectivity in a polymer molecule enforces an interdependence between features on different scales. As a consequence, the choice of where one building block ends and where the next one begins is not unique, and it is not obvious how to abstract from a fundamental degree of freedom and use it in an implicit way in a coarser model. Thus, we will use the generic term "coarse-grained" for any model employing the idea of soft interacting particles (blobs) equal to or larger than the monomers constituting the polymeric chains.
The degree of coarse-graining is application-driven and describes the number of atoms/molecules in a typical blob considered by the coarse-grained model. It is closely related to the minimal features of the atomistic model that should be retained in order to reproduce the desired properties from the coarse-grained model. Mapping an atomistic model to a coarse-grained one is very important in defining the positions of coarse-grained particles and directly influences the parameterization of the coarse-grained force field.
A general procedure in coarse-graining usually involves: defining the observable of interest and determining the degree of coarse-graining; deciding an appropriate mapping of the atomistic model to the coarse-grained one; deriving interactions between the coarse-grained particles; reproducing target functions with the coarse-grained model; optimizing parameters/functions in the coarsegrained model and validating its range of applicability; conducting coarse-grained simulations.

Mesoscopic Simulations
A major challenge in simulating realistic PNCs is that neither the length nor the time scales can be adequately addressed by atomistic simulations alone, because of the extensive computational load. Until relatively recently, a somewhat neglected level of description in materials modeling has been the mesoscopic regime, lying between atomic (or super-atomic like) particles and finite element-based representations of a continuum, and covering characteristic length scales of 10 −8 -10 −5 m. At this scale, the system is still too small to be regarded as a continuum, yet too large to be simulated efficiently using atomic models. In a more precise way, a mesoscale can be defined as an intermediate length scale at which the phenomena at the next level below (e.g. particle motions) can be regarded as having been equilibrated, and at which new phenomena emerge with their own characteristic time scales.
Among the several mesoscopic methods applied to the study of polymers, Self Consistent Field theory has been a well-founded tool [54]. This method adopts a field-theoretic description of the polymeric fluids and makes a saddlepoint (mean-field) approximation. An alternative to invoking the saddle-point approximation is performing a normal Metropolis Monte Carlo (MC) simulation, with the effective potential energy of the system given by field-theoretic functionals. One of the first attempts has been made by Laradji et al. [55] for polymer brushes and then by Daoulas and Müller [56] and Detcheverry et al. [57,58] for polymeric melts. The coordinates of all particles in the system are explicitly retained as degrees of freedom and evolve through MC moves. Tracking the motion of mesoscopic particles requires the use of stochastic dynamics [59].

Selected Unresolved Issues in PNCs
PNCs have been an area of intense industrial and academic research for the past twenty years. Irrespectively of the measure employed -articles, patents, or funding-efforts in PNCs have been exponentially growing worldwide over the last 10 years. PNCs represent a radical alternative to conventional filled polymers or polymer blends-a staple of the modern plastics industry. Considering the multitude of potential nanoparticles, polymeric resins, and applications, the field of PNCs is immense [60]. The restricted class of polymer nanocomposites defined above still presents a complex and fascinating problem in statistical mechanics due to the richness of physical phenomena in mixtures of flexible polymer coils and hard impenetrable objects. Despite the unprecedented efforts placed on PNCs research there are still open questions which have not been definitely addressed yet. In the following we will summarize a few of them; later we will analyze the perspective simulations and theoretical calculations have provided us with.
Fundamental issues and questions include, but are not limited to: the packing and structure of dense mixtures of long polymer chains and hard impenetrable fillers, in the presence of attractive, neutral or repulsive interactions; perturbation of polymer packing and the possible nonexistence of a bulk region of the polymer matrix, especially in the case of PNC films; non-universal filler-induced polymer conformational changes triggered by interfacial effects and/ or modification of the excluded volume screening mechanism of a pure polymer melt; the way in which geometric and chemical factors determine, in a nonadditive manner, the competing entropic and enthalpic contributions to the mixture free energy, miscibility and the physical nature of phase separated states. In all cases, the large particle surface-to-volume ratio leads to an amplification of a number of rather distinct molecular processes, implying pervasive interference between layers of polymers around nanoparticle surfaces.

Segmental Dynamics and the Glass Transition Temperature
When cooling a glass forming liquid, instead of freezing at a well defined temperature, one observes a huge increase of the viscosity which takes place continuously. Such glass formers can be either simple liquids or polymer liquids, and many features are similar in both regarding the glass transition. One defines the glass transition temperature, T g , as the temperature at which the dominant relaxation time on the molecular scale (or monomeric scale in the case of polymers) reaches about ~100 s, which corresponds typically to a viscosity of 10 12 Pa s in the case of simple liquids. Typically for such glass forming liquids, the viscosity increases by twelve orders of magnitude over a change of temperature of about 100 K down to T g . The underlying mechanisms involved in this dramatic increase are still poorly understood [61,62]. Experimental results on polymer dynamics and the glass transition in PNCs are not conclusive concerning the mechanism and the details of this modification. Increases or decreases in T g by as much as 30 K [63] have been reported depending on polymer-nanoparticle interactions. Reduction of T g has been reported in the case of weak interactions between filler and polymer [64]. In other cases the addition of nanoparticles causes no significant change to the glass transition of the polymer, presumably because effects causing increase and decrease of polymer mobility are present simultaneously, effectively canceling out each other [65]. Moreover, strong interactions between the filler particles and the polymer suppress crystallinity, yielding new segmental relaxation mechanisms in semicrystalline polymers, originating from polymer chains restricted between condensed crystal regions and the semi-bound polymer in an interfacial layer with strongly reduced mobility [66].
Concerning the spatial extent of the T g -shift, several studies [67,68] on PNCs show an increase of the glass transition temperature, suggesting that the mobility of the entire volume of the polymer is restricted by the presence of the nanoparticles. However, there are many experimental results suggesting that the restriction of chain mobility caused by the nanoparticles does not extend throughout the material but affects only the chains within a few nanometers of the filler surface [69]. The existence of such an interfacial layer seems relatively well-established in the case of silica-filled elastomers, however its exact nature is not well understood: experimental results have been described in terms of one or two distinct interfacial layers or a gradual change in dynamics with changing distance from the particle.

Enhancing Nanoparticle Dispersion by Surface Grafting
One of the biggest challenges is the rational control of filler clustering or aggregation, which often adversely affects material properties. The idea of achieving a good, uniform nanoparticle dispersion state has been the focus of considerable research, especially because of its favorable impact on optical and some mechanical properties of the resulting composites [70,71]. In the past few years, several research groups have modified the surface of nanoparticle fillers in an effort to improve their dispersion in a polymer matrix. A promising strategy for controlling the dispersion and morphology of PNCs is to graft polymer chains onto the nanoparticles to form a brush layer [33]. The free chain/brush interfacial interactions may be "tuned" by controlling grafting density, g , the degree of polymerization of the grafted chains, N g , and of the polymer host, N f , the nanoparticle size, p , and its shape. For example, if nanoparticles are grafted with chains compatible with the matrix polymer, filler dispersion is favored [72][73][74][75][76]. Motivated by this concept, experimentalists have synthesized nanometer sized particles with high surface grafting density [74,77,78]. At fixed polymer chemistry, when the molecular weight of matrix polymer is lower than that of grafted polymer, nanoparticles disperse. On the contrary, if the molecular weight of the matrix polymer is higher than that of the grafted polymer, nanoparticles are thought to aggregate [74]. Since both the matrix and the brush have the same chemical structure, the immiscibility for longer matrix chains is entropic in origin and attributable to the concept of "brush autophobicity" [72,[79][80][81][82][83]].

Mechanical and Rheological Properties
The dispersion of micro-or nano-scale rigid particles within a polymer matrix often-but by no means alwaysproduces an enhancement in the mechanical properties of these materials. As mentioned earlier, the most important application of this sort involves rigid inorganic particles (originally carbon black, later also silica) in a cross-linked elastomer matrix, where an improvement of mechanical properties is sought. This so-called rubber reinforcement is a complex phenomenon, which may involve an enhanced grip of tires on wet roads, an improved resistance to wear and abrasion, low rolling resistance, and an increase of tires' ultimate mechanical strength (toughness, tearing resistance).
There is a variety of phenomena seeking an explanation. For the sake of readability, we will focus on a subset of them. Under very small cyclic deformations, there is a linear viscoelastic regime characterized by a very significant increase (sometimes even by two orders of magnitude compared to the reference unfilled network) of the in-phase storage modulus, both under elongation and under shear [84]. At medium-to-large strains, filled elastomers exhibit a markedly non-linear response which is absent in unfilled elastomers ("Payne effect") [85]. The degree of non-linearity increases strongly with particle loading. An order-ofmagnitude drop in the modulus is often observed on going to 5-10 % deformation (under shear), bringing the asymptotic modulus of the filled systems much closer to that of the reference unfilled network.
Other related effects are commonly observed in filled elastomers. One is deformation hysteresis ("Mullins effect"): under cyclic deformation, the elastic modulus in the first cycle is higher than that in the following ones [86]. This points to some kind of "damage" of the material, which, however, is often reversible. The original properties can be recovered within a few hours, by high-temperature annealing of the sample. Secondly, fillers affect also the dissipative, out-of-phase components of the modulus. This is expected, since, probably, friction of the polymer chains against the filler surfaces, or of two particles against each other produces new energy dissipation mechanisms, which are absent in unfilled elastomers. Elastic and dissipative effects likely share a common origin. Finally, reinforcement effects have a remarkable temperature dependence. The small-strain (linear) modulus of filled rubbers decreases with temperature, pointing to important enthalpic effects. The situation is completely reversed compared to unfilled elastomers, where the modulus increases linearly with absolute temperature due to the entropic nature of rubber elasticity.

From Statistical Mechanics to Computer Simulations
Our discussion starts by introducing the formalism of statistical mechanics and briefly describing the basic methods used in computer simulations. We limit ourselves to the absolute minimum of definitions and methods to be presented, trying not to sacrifice consistency and rigor.

Motion in Phase Space
Statistical physics describes a system of N particles at a given state as one point in 6N/dimensional phase space, containing the atom positions and momenta (and neglecting the internal degrees of freedom of atoms) [87]. In classical mechanics, the state of the system is completely specified in terms of a set of generalized coordinates{ i } and generalized momenta{ i }, where i = 1, … , N [88]. We will refer to the 3N-dimensional set from which the generalized coordinates of the system { } ≡ 1 , 2 , … , N take on values as configuration space, and to the 3N-dimensional set from which the generalized momenta { } ≡ 1 , 2 , … , N take on values as momentum space. Any instantaneous microscopic state of the system can be written as a point: in the phase space of the system. The set of values of the macroscopic observables, such as temperature, pressure, etc., describes the system's macroscopic state. One macroscopic state corresponds to all the microscopic states that provide the same values of the macroscopic observables, defined by the macroscopic state. If we know the Hamiltonian, ℋ i , i , t , for the system, then the time evolution of the quantities i and i (i = 1, … , N) is given by Hamilton's equations of motion and where i = 1, 2, … , N and ∕ ≡ ∇ symbolizes the gradient operator with respect to the vectorial quantity . As the system evolves in time and its state changes, the system point traces out a trajectory in Γ-space. Since the subsequent motion of a classical system is uniquely determined by the initial conditions, it follows that no two trajectories in phase space can cross. If the Hamiltonian ℋ does not depend explicitly on time, then ℋ is a constant of the motion. Such is the case for conservative systems.

Time Average
Any property of the system, , is then a function of the points traversed by the system in phase space. The instantaneous property at a time t is (Γ(t)) and the macroscopically meaningful observable property obs is the time average of this, In experiments, the time average comes about quite naturally, since almost all experimental methods measure over much longer time scales than the longest relaxation time of the system. A straightforward approach, in order to get from molecular simulations, is to determine a time average, taking a discrete sum over M time steps of length Δt: This is the approach undertaken in Molecular Dynamics (MD) simulations, where the atoms' trajectory is followed as a function of time, so it is straightforward to obtain the average.

Phase Space Probability Density
When we deal with real systems, we can never specify exactly the state of the system, despite the deterministic character of its motion in phase space. There will always be some uncertainty in the initial conditions. Therefore, it is useful to consider Γ as a stochastic variable and to introduce a probability density (Γ, t) on the phase space. In doing so, we envision the phase space filled with a continuum (or fluid) of state points. If the fluid were composed of individual discrete points, then each point would be equipped with a probability in accordance with our initial knowledge of the system and would carry this probability for all time, since probability is conserved. Because state points must always lie somewhere in the phase space, we have the normalization condition where the integration takes place over the entire phase space. Similarly, the probability of finding the system in a small finite region D of Γ-space at time t is found by integrating the probability density over that region: The probability density for finding a system in the vicinity of Γ depends on the macroscopic state of the system, i.e. on the macroscopic constraints defining the system's size, spatial extent, and interactions with its environment. A set of microscopic states distributed in phase space according to a certain probability density is called an ensemble. A very important measure of the probability distribution of an equilibrium ensemble is the partition functionQ. This appears as a normalizing factor in the probability distribution defined by the ensemble.

Ensemble Average
The ergodic hypothesis, originally due to L. Boltzmann [89], states that, over long periods of time, the time spent by a system in some region of the phase space of microstates with the same energy is proportional to the volume of that region, i.e., that all accessible microstates are equiprobable over a long period of time. Ergodicity is based on the assumption (provable for some Hamiltonians) that any dynamical trajectory, given sufficient time, will visit all "representative" regions in phase space, the density distribution of points in phase space traversed by the trajectory converging to a stationary distribution.
According to the ergodic hypothesis we can calculate the observables of a system in equilibrium as averages over phase space with respect to the probability density of an equilibrium ensemble, ens (Γ). If ens (Γ) obeys the normalization condition, Eq. (6), on the entire phase space Γ and also is zero for all points outside the hypersurface ℋ(Γ) = E, the ensemble average can be defined as: In Monte Carlo (MC) simulations, the desired thermodynamic quantities are determined as ensemble averages: If we wish to obtain an average over points in phase space, there is no need to simulate any real time dependence of the system; one need only construct a sequence of states in phase space in the correct ensemble. In the context of equilibrium simulations, it is always important to make sure that the algorithm used in the simulation is ergodic. This means that no particular region in phase space should be excluded from sampling by the algorithm. Such an exclusion would render the simulation wrong, even if the simulated object itself is ergodic. From a practical point of view, the ergodicity of the system can and should be checked through reproducibility of the calculated thermodynamic properties (pressure, temperature, etc.) in runs with different initial conditions.

Microcanonical (NVE) Ensemble
In the microcanonical (NVE) ensemble the number of particles, N, the volume of the system, V and the total energy, E, are conserved. This corresponds to a completely closed system which does not interact in any way with the environment and lies in a container of fixed volume, V. For simplicity, we neglect the intramolecular degrees of freedom. Then, the system energy will be a sum of kinetic, , and potential, energies. Since the total energy E must be conserved, the criterion for adding states in the ensemble would be which means that not all, but only those states in phase space Γ that have total energy E 0 are allowed. This can also be stated so that the probability density of the ensemble is where is the Kronecker delta for a discrete system, and the Dirac delta function for a continuous one. The partition function in the microcanonical ensemble, Q NVE , is: The summation over states, ∑ Γ , is used if microscopic states are discrete and (Γ) has the meaning of a probability. For one-component classical systems, the sum can be replaced by an integral, yielding where N! takes care of the indistinguishability of particles of the same species and h 3N is the ultimate resolution for counting states allowed by the uncertainty principle.
The proper thermodynamic potential for the microcanonical ensemble is the entropy: where k B is the Boltzmann constant. We therefore have a statistical thermodynamic definition of entropy as a quantity proportional to the logarithm of the number of microscopic states under given N, V, E. Eq. (14) establishes a fundamental thermodynamic equation in the entropy representation.

Canonical (NVT) Ensemble
In the canonical (NVT) ensemble the number of particles, N, the volume of the system, V, and temperature, T are conserved. This corresponds to a closed system, which, however, can exchange heat with a large surrounding bath. The energy is fluctuating, but the temperature is constant, describing the probability distribution of energy fluctuations. The total energy of the system is given by its Hamiltonian, ℋ { i }, { i } . The probability density of the ensemble is: with k B being the Boltzmann's constant and Q NVT the partition function in the NVT ensemble: The thermodynamic function of the system is the Helmholtz energy: Eq. (17) defines a fundamental equation in the Helmholtz energy representation by expressing A as a function of N, V, T.

Isothermal -Isobaric Ensemble (NpT)
The isothermal-isobaric ensemble describes the equilibrium distribution in phase space of a system under constant number of particles, temperature, and pressure. The volume of the system is allowed to fluctuate. Thus, a point in phase space is defined by specifying V, i and i , where the domain from which the i s take on values depend on the value of V.
The probability density of the NpT ensemble can be derived from that of the microcanonical ensemble, by considering a bath around the system which acts as both a heat and a work reservoir for the system under study. The probability density, in a classical statistical mechanical formulation, is: where Q NpT is the isothermal-isobaric partition function: where V 0 denotes some basic unit of volume introduced to make the partition function dimensionless (the exact magnitude of V 0 is immaterial).
The connection between the formalism of the isothermal-isobaric ensemble and macroscopic thermodynamic properties is established via the Gibbs energy:

Configurational Integral
As long as the Born-Oppenheimer approximation [90] is valid (as it practically always is in equilibrium thermodynamics) the potential energy of the system, (Γ), depends only on the generalized coordinates, i . Similarly, the kinetic energy, (Γ), depends only on the momenta i . Hence we can rewrite the expression for the system Hamiltonian as: It can be now seen that, in a classical (as opposed to quantum mechanical) treatment, the partition function, e.g. of the NVT ensemble, factorizes into a product of kinetic (ideal gas) and potential (excess) parts: This can be written as a product of the ideal gas contribution, and the excess contribution as: where: is the so called configurational integral. The partition function of the ideal gas is: with Λ being the de Broglie or thermal wavelength: From the perspective of a particle-based model, the fundamental problem of equilibrium statistical mechanics, according to Chandler [91], is to evaluate a configurational partition function of the form of Eq. (24).
Two important consequences arise from Eq. (23). First, all the thermodynamic properties can be expressed as a sum of an ideal gas part and an excess part. The chemical details which govern the interactions between the atoms of the system are included in the latter. In fact, in MC simulations the momentum part of the phase space is usually omitted, and all calculations are performed in configuration space. The second important consequence of Eq. (23) is that the total average kinetic energy is a universal quantity, independent of the interactions in the system. Indeed, computing the average of with respect to the probability distribution of Eq. (15) and using the factorization of Eq. (23) we obtain that [87]: or, more generally ⟨ ⟩ = 1∕2N dof k B T for a system of N dof degrees of freedom. 1

Molecular Dynamics
In Cartesian coordinates, and under the assumption that the potential energy is independent of velocities and time, Hamilton's equations of motion read: hence where i is the force acting on atom i: with the gradient being taken keeping all positions other than i constant. Solving the equations of motion then involves the integration of 3N second-order differential equations Eq. (31) which are Newton's equations of motion.
The classical equations of motion possess some interesting properties, the most important one being the conservation law. If we assume that and do not depend explicitly on time, then it is straightforward to verify that Ḣ = dℋ∕dt is zero, i.e., the Hamiltonian is a constant of the motion. In actual calculations this conservation law is satisfied if there exist no explicitly time-or velocitydependent forces acting on the system. A second important property is that Hamilton's equations of motion are reversible in time. This means that, if we change the signs of all velocities, we will cause the molecules to retrace their trajectories backwards. The computer-generated trajectories should also possess this property. Concerning the solution of equations of motion, in the limit of very long times, it is clear that no algorithm provides an essentially exact solution. However, this turns out to be not a serious problem, because the main objective of an MD simulation is not to trace the exact configuration of a system after a long time, but rather to predict thermodynamic properties as time averages and calculate time correlation functions representative of the dynamics.
In the following we briefly describe the most popular family of algorithms used in MD simulations for the solution of classical equations of motion: the Verlet algorithms. Another family of algorithms comprises higher-order methods, whose basic idea is to use information about positions and their first, second, and higher order time derivatives at time t in order to estimate the positions and their derivatives at time t + Δt [93].
In general, higher-order methods are characterized by a much better accuracy than the Verlet algorithms, particularly at small times. However, their main drawback is that they are not reversible in time, which results in insufficient energy conservation, especially in very long-time MD simulations. On the contrary, the Verlet methods are not essentially exact for small times but their inherent time reversibility guarantees that the energy conservation law is satisfied even for very long times. This feature renders the Verlet methods, and particularly the velocity-Verlet algorithm, the most appropriate ones to use in long atomistic MD simulations.

Verlet Algorithm
The initial Verlet algorithm [94] ends up calculating the positions at time t + Δt by using two Taylor expansions around times t − Δt and t + Δt, respectively: Summing these two equations we obtain: The estimate of the new positions contains an error that is in the order of Δt 4 , where Δt is the time step employed in (33) our MD scheme. It should be noted that the Verlet algorithm does not use the velocities to compute the new positions. One can, however, derive the velocities from knowledge of the trajectory, using which is only accurate to order Δt 2 .

Velocity-Verlet Algorithm
The problem of defining the positions and velocities at the same time can be overcome by casting the Verlet algorithm in a different way. This is the Velocity-Verlet algorithm [95], according to which positions are obtained through the usual Taylor expansion whereas velocities are calculated through with all accelerations computed from the forces at the configuration corresponding to the considered time.

Langevin Dynamics
When a large system is simulated, it is generally desired to keep the number of degrees of freedom as low as possible. If a certain subset of particles can be distinguished, of which details of the motion are not relevant, these particles can be omitted from a detailed MD simulation. However, the forces they exert on the remaining particles must be represented as faithfully as possible. This means that correlations of such forces with positions and velocities of particle i must be incorporated in the equations of motion of particle i, while uncorrelated contributions can be represented by random forces. This brings us to the field of Langevin Dynamics (LD) [96,97]. In LD a frictional force, proportional to the velocity, is added to the conservative force, in order to mimic an implicitly treated background (e.g. solvent). The friction removes kinetic energy from the system. In order to compensate for the friction, a random force adds kinetic energy to the system. In the simplest case of LD, the random force is taken to have white-noise character, and no correlations between the various degrees of freedom are assumed to exist. Under these conditions, the velocity dependent frictional forces become proportional to the instantaneous velocity (36) of the particle considered. Thus, the equation of motion of a particle i is transformed into the stochastic equation: where the friction coefficient of a particle is denoted by i and the random force by i . The systematic force i is the explicit mutual force between the N particles of the system, which is to be derived from the potential (or free) energy of the system, which depends on the positions of all particles, denoted by i . The stochastic force, i (t), is assumed to be a stationary Gaussian random variable with zero mean and to have no correlations with prior velocities or with the systematic force: where ⟨… ⟩ ens denotes averaging over an equilibrium ensemble, indices , run over the Cartesian components (x, y and z), k B is Boltzmann's constant, T ref is the reference temperature of the LD simulation and W( ) is the (Gaussian) probability distribution of the stochastic force. According to van Gunsteren et al. [98], the solution of the linear, inhomogeneous, first order differential equation, Eq. (39), is:

Fluctuation-Dissipation Theorem
To generate a canonical ensemble, the friction and random force have to obey the fluctuation -dissipation theorem [99]. Einstein was the first to extract the diffusion coefficient and mobility in a special case of Brownian motion [100], and made allusions to the existence of a balance between random forces and friction. Then, Nyquist [101] formulated a limited version of the theorem, in his study of noise in resistors. Later, Callen and Welton [102] proved the theorem in a generalized form.
According to Kubo [103], two different kinds of the fluctuation-dissipation theorem can be distinguished. The fluctuation-dissipation theorem of the first kind relates the linear response of a system to an externally applied perturbation and a two-time correlation function of the system in the absence of external forces. The latter form is closely related to the famous Green-Kubo expressions for transport coefficients. The fluctuation-dissipation theorem of the second kind constitutes a relationship between the frictional and random forces in the system, relying on the assumption that the response of a system in thermodynamic equilibrium to a small applied perturbation is the same as its response to a spontaneous fluctuation [59].

Mori-Zwanzig Projection Operator Formalism
A formal way of deriving LD is the projection operator formalism of Zwanzig [104,105] and Mori [106,107]. The basis of the formulation is the assumption that we have partial knowledge of the evolving system, for example we can only track certain variables, while the effect of the other variables is modeled or approximated in a rigorous way. In this approach the phase space is divided into two parts, which we are called interesting and uninteresting degrees of freedom. For the approach to be useful, the uninteresting degrees of freedom should be rapidly varying in comparison to the interesting ones. Mori introduced two projection operators, which project the whole phase space onto the sets of interesting and uninteresting degrees of freedom. The full equations of motion are projected only onto the set of interesting degrees of freedom. The result is a differential equation with three force terms: a mean force between the interesting degrees of freedom, a dissipative or frictional force exerted by the uninteresting degrees of freedom onto the interesting ones and a third term containing forces not correlated with the interesting degrees. When the uncorrelated force is approximated by a random force the interesting degrees of freedom are considered independent of the uninteresting degrees of freedom.

Brownian Dynamics
If the friction exerted by the background to the particles under consideration is high, correlations in the velocity will decay in a time period over which changes in the systematic force are negligible. Such a system can be called overdamped. In this case, the left-hand side of Eq. (39) can be neglected, after averaging over short times. The result is Brownian Dynamics (BD), which is described by the position Langevin equation: The time scale separation makes possible the exchange of the second order stochastic differential equation Eq. (39) for a first order stochastic differential equation, Eq. (46), without affecting the dynamics on time scales longer than m i ∕ i .
Van Gunsteren and Berendsen [108,109] have proposed several algorithms for integrating Eq. (46). We will pay a closer look to the one which reduces to the Verlet algorithm for zero friction. If we assume a timestep of Δt, for large values of i ∕m i Δt in the diffusive regime, when the friction is so strong that the velocities relax within Δt, the BD algorithm reduces to: with i enumerating the particles, 1 ≤ i ≤ N, and marking a Cartesian component of the vectors. The components of the random displacement ℛ(Δt) are sampled from a Gaussian distribution with zero mean and width: The reader is reminded that the integration timestep Δt should be small enough, such that systematic forces do not change significantly over its duration. The integration scheme for BD Eq. (47) resembles a MC algorithm, except that there is no acceptance criterion. Rossky et al. [110] have derived the correct acceptance probability and introduced their method under the name "Smart Monte Carlo".

Dissipative Particle Dynamics
Molecular Dynamics (MD) is a powerful simulation technique capable of producing realistic results in a wide spectrum of applications. However, the computational cost of a detailed atomistic interaction model in that paradigm severely limits its applicability beyond extremely small spatial and short temporal scales. Within the family of simulation techniques designed to overcome the limitations of MD, we turn our attention to Dissipative Particle Dynamics (DPD), which allows the study of complex hydrodynamic phenomena in extensive scales. The DPD method was introduced in 1990s as a novel scheme for mesoscopic simulations of complex fluids [111,112]. In DPD simulations, the particles represent clusters of molecules that interact via conservative (non-dissipative), dissipative and fluctuating stochastic forces. Because the effective interactions between clusters of molecules are much softer than the interactions between individual molecules, much longer time steps can be taken relative to MD (47) simulations. This approach is ultimately based on the Langevin equation, the stochastic differential equation describing Brownian motion accounting for the omitted degrees of freedom by a viscous force and a noise term. The original DPD model tracks the equation of motion of the particles: where m i , i and i = i t are the mass, position and velocity of particle i, respectively. The total force, i , acting on each particle consists of three parts: where C ij , D ij and S ij represent the conservative, dissipative and stochastic forces between particles i and j, respectively. The conservative force depends on the distance between particles i and j, r ij and is directed along the unit vector of their separation, ̂ ij : where f C r ij is a non-negative (i.e. neutral or repulsive) scalar function determining the form of the conservative interactions, depending on the particular system of interest. In literature it is frequently implemented as a soft repulsion of the form: where ij is a parameter determining the maximum repulsion between the particles and r c is a cut-off distance.
The dissipative force, D ij , represents the effect of viscosity and depends on the relative positions and velocities of the particles. The form usually used for this interaction in DPD simulations is [113]: where is a friction coefficient, ij = i − j and w D r ij is a distance-dependent weighting function. The fluctuating random force depends on the relative positions of the particles, and is defined as: with being a coefficient, w S r ij is a distance-dependent weighting function and i is a random variable sampled from a Gaussian distribution with zero mean and unit variance. It should be noted that both the dissipative and the random force act along the particle separation vector and therefore conserve linear and angular momentum. Also, the resulting model fluids are Galilean invariant because the particle-particle interactions depend only on relative positions and velocities. The fluctuating stochastic force, S ij , heats up the system, whereas the dissipative force, D ij , reduces the relative velocity of the particles, thus removing kinetic energy and cooling down the system. Therefore, the stochastic and dissipative forces act together to maintain an essentially constant temperature which fluctuates around the nominal temperature of the simulation, T. Dissipative particle dynamics simulations can be thought of as thermostatted molecular dynamics simulations with soft particle-particle interactions.
Despite qualitative observations, there was no theoretical justification that DPD produces the correct hydrodynamic behavior until Español and Warren [114] formulated the Fokker-Planck equation for studying the equilibrium properties of the stochastic differential equation describing DPD. Later, Español [115] derived the macroscopic hydrodynamic variables starting from the microscopic description. In order to recover the proper thermodynamic equilibrium for a DPD fluid at a temperature T, the coefficients and the weighting functions for the dissipative and random forces should be related by: and as required by the fluctuation-dissipation theorem. All interaction energies are expressed in units of k B T, which is usually assigned a value of unity. One straightforward and commonly used choice is: where r c is the cut-off distance of the the dissipative and the random force. In conventional DPD formulation, it usually takes the same value as the cut-off distance of the conservative force but it can vary in order to modify the dynamic properties in DPD simulations. For conventional DPD simulations, the exponent of the weighting function, s, is set equal to 2 with w D and its gradient being continuous at r ij ∕r c = 1.
Summarizing, Español and Warren [114] established a sound theoretical basis for DPD and Groot and Warren [116] obtained parameter ranges to achieve a satisfactory compromise between speed, stability, rate of temperature equilibration, and compressibility. Unlike traditional DPD methods using conservative pairwise forces between particles, the multi-body DPD model presented by Pagonabarraga and Frenkel [117] assumes that the conservative force depends on the instantaneous local particle density, which in turn depends on the positions of many neighboring particles. As far as the integration of the DPD equations of motion is concerned, Pagonabarraga et al. [118] proposed a leap-frog scheme which was self-consistent and satisfied a form of microscopic reversibility. Thus, the correct equilibrium properties could be recovered from trajectories generated with that algorithm.

Monte Carlo
The Monte Carlo (MC) method is a statistical approach for finding approximate solutions to problems by means of random sampling. In addition to molecular simulations and physics, it is widely applied in other natural sciences, mathematics, engineering and social sciences [119]. The earliest treatments in the subject, such as this by Babier [120], were made in connection with the "Buffon's needle problem". 2 According to Metropolis [121], the invention of the modern class of MC algorithms is due to Enrico Fermi, when he was studying the properties of the then newly-discovered neutron in 1930. It was further developed during the 1940s by physicists working in the nuclear weapons program of the United States, at the Los Alamos National Laboratory. The technique was given its name by Nicholas Metropolis, in reference to the famous casino in Monaco, considering the use of randomness and the repetitive nature of the sampling process.
In their simplest version, MC simulations of simple fluids are carried out by sampling trial moves for the molecules from a uniform distribution. For example, in a canonical (NVT) ensemble simulation, a molecule is chosen at random, and then displaced, also randomly to a new position. The trial move is accepted or rejected according to an importance sampling scheme [93,122,123]. A frequently used importance sampling algorithm is the Metropolis algorithm, originally derived for the specific case of the Boltzmann distribution [122] and later generalized to other distributions [124] which need not to be analytical (e.g. the force-bias method of Pangali et al. [123] provides a classical example of such algorithms).
The probability of accepting a move, P accept , of the form: will asymptotically sample the configuration space according to a probability P. In Eq. (58), P accept is the probability with which trial moves are accepted or rejected, P( | ) is the transition probability of making a trial move from state to state , and P( ) is the probability of being at state . This means that, at equilibrium, the average number of accepted trial moves that result in the system leaving state must be exactly equal to the number of accepted trial moves from all other states to the state . This is a looser statement of the detailed balance condition, reflected in Eq. (58), that at equilibrium the average number of accepted moves from to any other state should be exactly canceled by the number of reverse moves.
In the original Metropolis scheme [122], the probabilities P( | ) form a symmetric matrix, constructing a Markov chain that has the Boltzmann distribution as its equilibrium distribution. In this case, there is no bias involved in making the move and Eq. (58) reduces to the standard Metropolis acceptance criterion: The advanced MC methods are based on judicious choices of P( | ) [93]. It should be noted that the simulation steps in the MC technique are steps in configuration space and there is no notion of "time" in MC simulations. This is contrast to MD, where the simulation steps are explicit time steps. Moreover, a computational advantage of MC over MD is that only the energy needs to be calculated, not the forces, rendering the Central Processing Unit (CPU) time needed per step smaller than that of an MD simulation.

Reduced Units
Molecular simulations can conveniently be performed in non-dimensional or reduced units, based on the characteristic physical dimensions of the system under study. Working with reduced units is preferred mainly because they are physically easier to interpret, and the results obtained become applicable to all materials modeled by the same potential. Reduced units are obtained by expressing all quantities of the simulation in terms of selected base units which are characterizing the system, in order to make equations dimensionless. Table 1 presents some reduced quantities. For example, in the case of the Lennard-Jones potential, the particle diameter, , the depth of the potential well, , together with the mass of the simulated particles, m, provide a meaningful and complete set of base units for simulations.

Chain Dimensions in the Bulk
One of the most important and probably the most fundamental question in the area of PNCs is how the size of the polymer chains is affected by the dispersion of nanoparticles. There has been considerable controversy in the experimental literature over whether nanoparticles cause chain expansion (swelling) [125,126], contraction, [127] or neither [128][129][130][131][132][133]. The sign (attractive or repulsive) and strength of the nanoparticle/polymer interactions, the relative dimensions of the chains with respect to the size of the nanoparticle, R g ∕R n , and the exact state of dispersion, have been identified as the key factors that can account for the aforementioned differences in the structure of the matrix chains.

Experimental Findings
Chain conformations in PNCs have been primarily measured by small angle neutron scattering (SANS). These measurements are greatly facilitated by combining deuterated and hydrogenated chains such that the average scattering length density of the polymer blend matches that of the filler. This zero average contrast condition [134], which is hard to achieve, minimizes the scattering due to the nanoparticles. To date, studies which report increases in polymer dimensions, in the case of spherical nanoparticles, invoke the presence of attractive nanoparticle/polymer interactions, combined with R g > R n , and good nanoparticle dispersion [135], to conclude that the nanoparticles behave as a good solvent for the polymer chains. However, even though the existence of a shell containing polymer of reduced mobility is acceptable in nanocomposites composed of strongly interacting particles and polymer, e.g. composed by silica and PMMA, the size of the chains, e.g. in terms of their radius of gyration, R g , is intrinsically independent of the the volume fraction, , (up to ≃ 0.1) and the polymer-to-particle size ratio [132]. All other studies on spherical nanoparticles showed little if any changes in polymer R g , that is where the nanoparticle-polymer intractions are believed to be athermal, or significant nanoparticle aggregation was present, due to unfavorable nanoparticle/polymer interactions [136].
In an early study of a poly(dimethylosiloxane) / polysilicate (R n = 1 nm) nanocomposite [125], a significant increase of the polymer chain dimensions (reaching 60% expansion at nanoparticle volume fraction, = 40%) was observed for R g > R n and a decrease in polymer dimensions for R g ≤ R n . Neutron scattering studies of an athermal polystyrene (PS) PNC, indicated swelling of the matrix chains, induced by dispersed tightly cross-linked PS nanoparticles [126]. PS chains around crosslinked PS particles (R n = 2 − 3.6 nm) were found to be expanded by up to 20 % relative to their unperturbed size, when their unperturbed radius of gyration was comparable to or larger than the radius of the dispersed particles. More recent studies of PS/ silica nanocomposite [128,129,131] for R g ∕R n = 2 − 4, [131] and poly(ethylene-propylene)/silica nanocomposites, [127] found no perturbation of the matrix chain dimensions.
We may attribute the qualitatively different trends deduced by different experimental studies to several factors, including, but not limited to the following: (a) most of the polymers used exhibit significant polydispersity, (b) particle dispersion/agglomeration cannot be adequately quantified, (c) the molecular weight of the isotopic polymers blended with the filler was quite different in at least one case. The compound effect of these factors can result in significant uncertainty in the chain dimensions measured. Molecular simulations can shed some light on the role of nanoparticles on chain dimensions, especially in regimes where it is hard to conduct experiments.

Insight Obtained from Simulations
From the point of view of molecular simulations, there also exists controversy as to whether the incorporation of nanoparticles in a polymer melt causes polymer chains to expand, remain unaltered or reduce [137][138][139] their dimensions compared to their size in the bulk material. To date, all studies have indicated that, irrespectively of the absolute dimensions of the chains in the interparticle region, these retain their unperturbed Gaussian scaling. This is a striking feature, resembling the scaling behavior of chains in thin films [140,141], where chain conformations parallel to the surface assume their unperturbed values even for film thicknesses < R g . Most of the simulation works have addressed polymer dimensions in nanocomposites below the percolation threshold ( c = 31% [142]), except the early works of Vacatello [137][138][139]143] that were implemented at constant density and spatially frozen nanoparticles.
Sen et al. [128] employed polymer reference interaction site model (PRISM) [144,145] calculations in order to interpret small angle neutron scattering findings on polystyrene loaded with spherical silica nanoparticles under contrast-matched conditions. They considered blends with 66 wt% hPS and 34% dPS, which almost contrast match the silica. Nanocomposites with 0, 2.9, 6.1 and 9.1% vol silica were prepared for each molecular weight and 15.9 and 27.4 vol% for the higher molecular weights considered. However, in their experiments, as in earlier studies [146,147], the particles were imperfectly mixed with the polymeric matrix, with particles being surrounded by "voids", especially at large filler contents. In parallel, the PRISM theory was applied, by modeling the fillers as hard spheres and polymers as freely jointed chains with a realistic persistence length. Polymer-polymer and particle-particle interactions were taken to be hard-core, while monomers and filler interact via an exponentially decreasing attraction over a predefined spatial range. From the experimental point of view (Fig. 1), the low-q intensity increases dramatically with increasing silica content, especially for loadings ≤10 vol%, implying that the matrix is not totally contrast matched to the filler (unsurprising in light of voids surrounding particles [146]). However, the scaling and dimensions of the polymer chains can be obtained from the high-q intensity which is expected to be independent of the filler structure [146]. Their results (Fig. 1) showed that chain conformations follow unperturbed Gaussian statistics independent of chain molecular weight and filler composition. Liquid state theory calculations were consistent with this conclusion and also predicted filler-induced modification of interchain polymer correlations which had a distinctive scattering signature that was in nearly quantitative agreement with the experimental observations. The chain R g varied from ~8 nm (90 kg/mol) to 22 nm (620 kg/mol), bracketing the nanoparticle diameter (~14 nm), suggesting that the ratio of the particle size to R g was not an important variable in that context.
The structure of a polystyrene matrix filled with tightly cross-linked polystyrene nanoparticles, forming an athermal nanocomposite system, has been investigated by means of a Monte Carlo sampling formalism by Vogiatzis et al. [148]. Although the level of description is coarsegrained (e.g., employing freely jointed chains to represent the matrix), the approach developed aims at predicting the behavior of a nanocomposite with specific chemistry quantitatively, in contrast to previous coarse-grained simulations. A main characteristic of the method was that it treats polymer-polymer and polymer-particle interactions in a different manner: the former are accounted for through a suitable functional of the local polymer density, while the latter are described directly by an explicit interaction potential. The simulation methodology was parameterized in a bottom-up fashion in order to mimic the experimental studies. Many particle systems, with volume fractions up to 15 vol%, were simulated. The positions of the nanoparticles were held constant in the course of the simulation, while polymeric chains were allowed to equilibrate via a combination of MC moves. The generation of many independent initial configurations compensated for the immobility of the particles along the simulation. The values of the radius of gyration R g , relative to the value for the pure polymer melt R g0 , are shown in Fig. 2 as a function of the nanoparticle volume fraction for the four different chain molecular weights used in that work (23, 47, 93 and 187 kg/mol). In general, an expansion of polymeric chains with increasing nanoparticle volume fraction can be observed for all chain lengths. This expansion is maximal for 23 kg/mol, where the unperturbed radius of gyration R g,0 ≃ 42 Å is comparable to the radius of the nanoparticle, R n = 36 Å. It seems that there is a tendency of chains to swell when their dimension is equal to or approaches the dimension of the nanoparticle. This observation is in very good quantitative agreement with experimental data reported for the same system [126]. In all other cases, the swelling due to the presence of the nanoparticles could hardly be discerned.
Karatrantos et al. [135] have investigated the effect of various spherical nanoparticles on chain dimensions in polymer melts for high nanoparticle loading which was larger than the percolation threshold, using coarse-grained molecular dynamics simulations of the Kremer-Grest model [149]. Their results, presented in Fig. 3, revealed different behavior of the polymer chains in the presence of repulsive or attractive particles. In nanocomposites containing repulsive nanoparticles (black symbols), the polymer dimensions were not altered by the particle loading. These authors reported that the polymers were phase separated from the repulsive nanoparticles (of R n = 2) in the nanocomposites, thus, there were no changes in the radius of gyration values. On the contrary, in the nanocomposites containing attractive nanoparticles, the overall polymer dimensions increased dramatically at high particle loading. In particular, the magnitude of expansion of polymer dimensions was larger for polymers with N = 200 following qualitatively the experimental data [125,150]. The relation R g ∕R g,0 = (1 − ) −1∕3 , included in Fig. 3 was proposed by Frischknecht et al. [151] for predicting the polymer expansion due to the excluded volume introduced by the nanoparticles, assuming no change in density on mixing. Finally, Karatrantos et al. [135] reported that polymer chains, in all cases considered, did not depart from Gaussian statistics.
Mathioudakis et al. [153] applied a hierarchical simulation approach in order to study the behavior of PS-SiO 2 nanocomposites. Two interconnected levels of representation were employed. (a) A coarse-grained one [154], wherein each polystyrene repeat unit was mapped into a single "superatom" and each silica nanoparticle into a sphere. The smoother effective potential energy hypersurface of the coarse-grained representation permitted its equilibration at all length scales by using powerful connectivity-altering Monte Carlo algorithms [155]. (b) A united atom representation, wherein polymer chains were represented by a united-atom model and a silica nanoparticle was represented in full atomistic detail. Coarse-graining and reverse-mapping between the two levels of representation was accomplished in a manner that preserved tacticity and respected the detailed conformational distribution of chains [156]. At the coarser level, these authors estimated the root-mean-square radius of gyration as a function of the molecular weight for neat and nanocomposite polystyrene systems. Their results are presented in Fig. 4 along with neutron scattering results for bulk monodisperse PS [152] from 21 to 1100 kg/mol. As far as the nanocomposite polystyrene systems are concerned, the presence of the nanoparticles affected the root-meansquare radius of gyration only slightly.

Experimental Findings
SANS measurements show a clear scattering signature of a polymer bound layer around the particles, which arises due to a scattering length density different from the bulk polymer matrix material, either due to H or D enrichment or a modification of the polymer density in the bound layer compared to the surrounding polymer matrix [132,133]. The measurements of Jouault et al. [133] revealed that the bound layer is independent of the particles' volume fraction. Then, as observed by Jiang et al. [157], the bound layer volume fraction is larger at the surface (that region being mostly composed of loops) and decreases at larger distances as the bound layer becomes more diffuse due to the contribution from the tails. One can estimate the thickness of the bound polymer layer around 2 nm. However, this thickness value is a simplification because it does not completely describe the complex chain behavior, some aspects of which will be analyzed below.

Insight Obtained from Simulations
The local density of the polymer in the proximity of the surface of a filler is often employed as an indication of the strength of polymer-surface interactions and a decrease of the first peak of the radial density distribution is expected with curvature [158]. At this point we resort to the detailed analysis of Pandey and Doxastakis [159] concerning a polyethylene layer close to a filler surface (Fig. 5). These authors coupled the application of preferential sampling techniques [160] with connectivity-altering Monte Carlo algorithms [161,162] in order to explore the configurational characteristics of a polyethylene melt in proximity to a silica surface or around a nanoparticles and the changes induced by high curvature when the particle radius is comparable to the polymer Kuhn segment length. The inset to Fig. 5 shows that indeed as we move from a flat surface to a smaller nanoparticle a decrease is observed with the exception of the fullerene where a significantly higher density is found. To investigate further the concentration of adsorbed monomers, these authors followed the use of a simple distance criterion (adsorbed polymer chains have an atom within 0.6 nm of an atom of the surface; introduced by Daoulas et al. [163]) to decompose polymer segments according to the Scheutjens-Fleer theory [164] into trains, tails and loops. Tails are the segments which are hinged to the surface at one end while the other end is dangling freely into the bulk polymer. Train segments consist of monomers consecutively adsorbed on the surface. The loop segments are constituted by the monomers in-between two train segments, which are not adsorbed on the surface. Figure 5 exhibits three regimes for adsorbed chains: a first layer of adsorbed monomers constituting train segments, a second layer where a decay is dominated by a decrease of loop segments while tail density is constant and a third layer where tail segments extend in the bulk melt. As shown in the inset to Fig. 5(a) the area under the first peak broadens as we move to smaller particles.
An interesting feature of the interfacial systems to study is the number of monomers that are in contact with the surface. To this end, Pandey and Doxastakis [159], defined a surface concentration by integrating the density profile of train segments: where R n is the radius of the nanoparticle, train is the density profile of train segments, and A is the accessible surface area to the polymer. If we assume that nanoparticles are spheres surrounded by a constant density of polymer, 0 , in a layer of Δr thickness, the surface concentration is given by: where a constant density is multiplied with the ratio of the volume of a spherical shell representing the first adsorbed monolayer to the surface of the sphere. The geometric predictions following the above line of reasoning, are shown for different chain lengths by the dashed lines in Fig. 5(b). It is apparent that a modest increase and ultimate leveling 5 a Density distribution of a polyethylene melt as a function of distance from the surface of a filler (graphite slab, silica nanoparticle or fullerene C 60 ). The decomposition into tails, trains and loops is carried out following Scheutjens and Fleer [164]. The inset provides profiles for selected systems. b Surface concentration together with predictions based upon geometrical arguments for ideal spheres and surface monomer density in the proximity of silica slabs. (Reprinted from [159], with the permission of AIP Publishing.) off of the surface concentration with decreasing nanoparticle radius is observed in sharp contrast to the estimations based on the geometric arguments, which predict a continuous increase. The lower than anticipated increase of surface concentration with curvature suggests that collective properties beyond the enthalpic interactions appear to play a crucial role on surface concentration. At the extreme limit where particles are comparable to the polymer Kuhn segment length, curvature penalizes the formation of long train segments. As a result, an increased number of shorter contacts belonging to different chains were made, competing with the anticipated decrease of the bound layer thickness with particle size if polymer adsorbed per unit area remained constant. Starr et al. [165] carried out extensive molecular dynamics simulations of a single nanoscopic filler particle surrounded by a dense polymer melt. The polymers were modeled as chains of monomers connected via a finitely extensible nonlinear elastic (FENE) anharmonic potential and interacting via a Lennard-Jones potential. That type of "coarse-grained" model is frequently used to study general trends of polymer systems but does not provide information for a specific polymer. The filler particle shape was icosahedral with interaction sites assigned at the vertices, at four equidistant sites along each edge, and at six symmetric sites on the interior of each face of the icosahedron. These authors considered a filler particle with an excluded volume interaction only, as well as one with excluded volume plus attractive interactions in the dilute nanoparticle regime (where bulk chain dimensions are unlikely to be affected by the confinement between nanoparticles). By focusing on the dependence of R g on the distance d from the filler surface, these authors were among the first to report a change in the overall polymer structure near the surface. In Fig. 6, R 2 g , as well as the radial component of from the filler center R ⟂2 g (approximately the component perpendicular to the filler surface) for both attractive and nonattractive polymer-filler interactions at one temperature. A striking feature of Fig. 6 is that R 2 g increases by about 30% on approaching the filler surface, while at the same time R ⟂2 g decreases by more than a factor of 2 for both (attractive and neutral) systems. The combination of these results indicates that the chains become slightly elongated near the surface and flatten significantly. The range of the flattening effect roughly spans a distance of an unperturbed radius of gyration, R g,0 , from the surface and depends only weakly on the simulation temperature, T. Moreover, the chains retain a Gaussian structure near the surface.
Mathioudakis et al. [153] studied the shape of chains in the presence of a silica nanoparticle by employing coarsegrained MC simulations. These authors analyzed the eigenvalues of the the radius of gyration tensor, which serve as a measure for characterizing the shape of the polymer chains. In the polymer melt, the intrinsic shape of chains is that of a flattened ellipsoid or soap bar [166]. Following ref. [166], Mathioudakis et al. diagonalized the instantaneous radius of gyration tensor of every chain to determine the eigenvalues L 2 3 ≥ L 2 2 ≥ L 2 1 (squared lengths of the principal semiaxes of the ellipsoid representing the segment cloud of the chain) and the corresponding eigenvectors (directions of the principal semiaxes). The three semiaxes are generally unequal. The sum L 2 1 + L 2 2 + L 2 3 equals the squared radius of gyration of the chain. These authors observed that, when the distance of the center of mass of the chain from the center of the nanoparticle was shorter than the mean size of the chain, the chains expanded along their principal semiaxis, L 3 . That led to an increase of the radius of gyration, R 2 g = L 2 1 + L 2 2 + L 2 3 near the nanoparticle. The deformation of the molecules was smaller for chains whose dimensions exceed by far the radius of the nanoparticle. Far from the surface of the nanoparticle, the sum of the squares of the principal semiaxes (sum of the eigenvalues of the radius of gyration tensor) reaches the bulk average value of the squared radius of gyration of PS, since the molecules were not affected by the presence of the nanoparticle. These results are shown in Fig. 7. Changes in the The increase of R g , coupled with the decrease of R ⟂ g , indicates that the chains become increasingly elongated and flattened as the surface of the particle is approached. (Reprinted figure with permission from [165]. Copyright 2001 by the American Physical Society) intrinsic shape of chains were quantified as a function of distance of the center of mass of the chain from the center of mass of the silica particle by calculating the ratio of largest to smallest eigenvalue of the radius of gyration tensor. This local anisotropy of the chains as a function of distance from the centre of mass of the nanoparticle is also shown in Fig. 7.
The presence of the filler surface also influences the orientation of the chains. Ndoro et al.
[158] studied the distance dependence of the angle between the longest axis of the radius of gyration tensor and the surface normal of bare silica nanoparticle. Their results are presented in Fig. 8. The observation that the free polymer chains generally prefer to align tangentially to the ungrafted surface is in agreement with conclusions from other researchers [148,167]. In their coarse-grained model using Monte Carlo simulations, Vogiatzis et al. [148] studied the orientational angles of local chain segments. They also concluded that chain segments in the vicinity of the nanoparticle surface were structured and oriented tangentially to the interface.
Bačová et al. [168] conducted atomistic molecular dynamics simulations of graphene-based polymer nanocomposites composed of hydrogenated and carboxylated graphene sheets dispersed in polar and nonpolar short polymer matrices, in order to gain insight into the effects of the edge group functionalization of graphene sheets on the properties of hybrid graphene-based materials. Poly(ethylene oxide) and polyethylene serve as the polar and nonpolar matrix, respectively. In Fig. 9 the structural properties of the short polymer chains, i.e. their mean square end-to-end distance ⟨ R 2 e ⟩ , and the radius of gyration, , for the chains, whose center of mass is placed within a given layer. The five layers employed are set up in accordance with the positions of the minima in the density profiles (cf. Fig. 9). The results for all five parallel layers and both polymer matrices are plotted in Fig. 9. The data are normalized by the bulk values. The error bars correspond to the standard deviation and were obtained through

End Grafted Polymers onto Nanoparticles
Controlling the spatial dispersion of nanoparticles is critical to the ultimate goal of producing polymer nanocomposites with desired macroscale properties. Experimental studies [14,73,150] have shown that, often, nanoparticles tend to aggregate into clusters, with the property improvements connected to their nanoscale dimension being lost. One commonly used technique for controllably dispersing them is end grafting polymer chains to the nanoparticle surface, so that nanoparticles become "brush coated" [14]. When the coverage is high enough, the nanoparticles are sterically stabilized, which results in good spatial dispersion [173,174]. Moreover, spherical nanoparticles uniformly grafted with macromolecules robustly self-assemble into a variety of anisotropic superstructures when they are dispersed in the corresponding homopolymer matrix [14].
A simpler system that is useful for understanding polymer brushes grafted on spherical nanoparticles immersed in melts is that of a brush grafted to a planar surface in Fig. 9 Mean square end-to-end distance of polymer chains located at different layers parallel to graphene, normalized by the same quantity measured in bulk. In the inset to the figure, the normalized radius of gyration is plotted for the same set of data. The first layer extends within distances 0.0 and 0.6 nm from the graphene sheet, the second between 0.6 and 1.0 nm, the third between 1.0 and 1.5 nm, the fourth between 1.5 and 2.0 nm and the last one, fifth, between 2.0 and approximately half the edge length of the simulation box, 5.0 nm. (Reprinted with permission from Ref. [168]. Copyright (2011) American Chemical Society.) contact with a melt of chemically identical chains. Important molecular parameters for this system are the Kuhn segment length of the chains, b, the lengths (in Kuhn segments) of the grafted, N g and free, N f , chains, and the surface grafting density (chains per unit area), . The case of planar polymer brushes exposed to low molecular weight solvent was studied theoretically by de Gennes [175] and Alexander [176]. These authors used a scaling approach in which a constant density was assumed throughout the brush: all the brush chains were assumed to be equally stretched to a distance from the substrate equal to the thickness of the brush. Aubouy et al. [177] extracted the phase diagram of a planar brush exposed to a high molecular weight chemically identical matrix. Their scaling analysis is based on the assumption of a steplike concentration profile and on imposing the condition that all chain ends lie at the same distance from the planar surface. Five regions with different scaling laws for the height, h, of the brush were identified. For low enough grafting densities, < N −1 g a −2 (with a being the monomer size) and short free chains, N f < N 1∕2 g , the brush behaves as a swollen mushroom, with h ∼ N 3∕5 g . By keeping the grafting density below N −1 g a −2 and increasing N f , so that N f > N forces the chains to stretch, thus leading to a brush height scaling as h ∼ N g .
Wijmans and Zhulina [178] employed similar ideas in order to understand the configurations of polymer brushes grafted to spherical nanoparticles dispersed in a polymer melt. Here the radius of the nanoparticle, R n , enters as an additional parameter. Long polymers grafted to a surface at fixed grafting density, , are strongly perturbed from their ideal random-walk conformation [179]. Planar geometry scaling (infinite radius of curvature) is inadequate to explain the case when the particle size is reduced to a level comparable with the size of the brushes. The SCF theory has been applied to convex (cylindrical and spherical) surfaces by Ball et al. [179]. For the cylindrical case, under melt conditions, it was found that the free ends of grafted chains are excluded from a zone near the grafting surface. The thickness of this dead zone varies between zero for a flat surface to a finite fraction of the brush height, h, in the limit of strong curvature, when R n ∕h is of order unity, with R n being the radius of curvature of the surface.
Borukhov and Leibler [180] presented a phase diagram for brushes grafted to spherical particles, in which the five regions of the work of Aubouy et al. can still be located, but they also provide the scaling of the exclusion zone, where matrix chains are not present. Such exclusion zones have been observed in a special limiting case of grafted polymers, namely, star shaped polymers. Daoud and Cotton [181] showed that, in the poor-solvent limit, the free ends of the chains are pushed outward, because of high densities near the center of the star. The Daoud-Cotton model assumes that all chain ends are a uniform distance away, while the Wijmans-Zhulina model [178] has a well-defined exclusion zone. For Θ solvents, in the limit of large curvature (small particle radius, R n ), the segment density profile, (r), decreases with the radius as [178] (r) ∝ 1∕2 R n ∕r . It must be noted that is not linear in because the brush height depends on . In the limit of small curvature (large R n ), a distribution of chain ends must be accounted for [182], leading the segment density profile to a parabolic form: [178] (r � ) = where b is the statistical segment length, r � = r − R n is the radial distance from the particle surface, h 0 is the effective brush height for a flat surface and h is the brush thickness. For large nanoparticles, the above form asymptotically recovers the planar result (h → h 0 ). In the case of intermediate particle radii, a combination of large and small curvature behaviors is anticipated: [178] the segment density profile exhibits large curvature behavior near the surface of the particle, followed by a small curvature behavior away from it. Finally, following Daoud and Cotton, the brush height is expected to scale as h ∝ 1∕4 N 1∕2 g . Recently, Chen et al. [183] revisited the scaling laws for spherical polymer brushes and identified significant assumptions overlooked by Daoud and Cotton.

Experimental Findings
Hasegawa et al. [173] used rheological measurements and SCF calculations to show that particles are dispersed optimally when chains from the melt interpenetrate, or wet, a grafted polymer brush ("complete wetting"). This occurs at a critical grafting density, which coincides with the formation of a stretched polymer brush on the particle surface [175,176]. Grafting just below this critical density produces aggregates of particles due to attractive van der Waals forces between them. The results of these authors suggest that mushrooms of nonoverlapping grafted polymer chains have no ability to stabilize the particles against aggregation ("allophobic dewetting"). Grafting just above this critical density also results in suboptimal dispersions, the aggregation of the particles now being induced by an attraction between the grafted brushes. For high curvature (small radius) nanoparticles, the polymer brush chains can explore more space, resulting in less entropic loss to penetrate the brush, reducing the tendency for autophobic dewetting.
Until recently, the experimental verification of theoretical and simulation predictions was mostly limited to global information concerning the brush, such as its average height, but not its profile [184]. Recently Chevigny et al. [185] used a combination of X-ray and Small Angle Neutron Scattering (SANS) techniques to measure the conformation of chains in polystyrene polymer brushes grafted to silica nanoparticles with an average radius of 13 nm dispersed in polystyrene matrix. They found that, if the molecular weight of the melt chains becomes large enough, the polymer brushes are compressed by a factor of two in thickness compared to their stretched conformation in solution. Also, polymer brushes exposed to a high molecular weight matrix are slightly compressed in comparison to brushes exposed to a low molecular weight matrix environment. This observation implies a wet to dry conformational transition. The low molar mass free chains can penetrate into the grafted brush and swell it ("wet" brush). Conversely, when grafted and free chain molar masses are comparable, free entities are expelled from the corona ("dry" brush). Later, they examined the dispersion of these grafted particles in melts of different molar masses, M f [186]. They showed that for M g ∕M f < 0.24, the nanoparticles formed a series of compact aggregates, whereas for M g ∕M f > 0.24, the nanoparticles were dispersed within the polymer host.

Insight Obtained from Simulations
Klos and Pakula [187] simulated linear flexible polymers of N g repeat units, end-grafted at density onto a spherical surface of radius R ("hairy nanoparticle"), including the case of flat impenetrable wall (R → ∞) using their cooperative motion algorithm [188,189]. The simulations were carried out for a wide range of parameters characterizing the hairy surfaces (N g , , and R) and concerned in detail the influence of length of matrix chains on the anchored ones. That was achieved by gradually varying the polymerization degree N of the matrix chains between the two extremes of a dense melt of identical chains (N = N g ) and a simple solvent consisting of single beads (N = 1). Their analysis of free grafted chain-ends concentrations, fe (z), is shown in Fig. 11 (a, b) for R =3, 10 and for =1, 0.2, respectively. The length unit used in their work was c / 2 with c being the lattice constant of the employed lattice Monte Carlo simulations. The curves indicate how the medium in which the hairy sphere is immersed influences the profiles of free ends of the grafted chains. For both sizes of the spheres, the observed tendency is such that the longer the matrix chains become, the closer to the surface the free ends concentrate. In particular, this is also the case for chains grafted to a flat surface, as presented in Fig. 11(c). Furthermore, for N = N g , R = 10, and R → ∞, the concentration of the free ends is finite even at the surface, which means that a small fraction of the ends concentrate in that region, creating grafted chain loops, in agreement with earlier Molecular Dynamics simulations of brushes on flat surfaces by Grest [190].
Spatial integration of the radial mass density profiles around the nanoparticle allows for estimating the height of the grafted polymer brush, which is usually defined as the second moment of the segment density distribution, (r), as [178,191]: with respect to the height h ≡ r − R n . However, comparison with experimental brush heights requires a measurement of where the major part of the grafted material is found. To this effect, the brush height can also be arbitrarily defined as the radius marking the location of a spherical Gibbs dividing surface, in which 99% of the grafted material is included. The theory of spherical polymer brushes was pioneered by Daoud and Cotton [181]. In analogy to the scaling model developed by Alexander and de Gennes for planar interfaces, Daoud and Cotton developed a model for spherical surfaces through geometric considerations based on starlike polymers. The spherical brush is divided into three regions, an inner meltlike core region, an intermediate concentrated region (dense brush), and an outer semidilute region (swollen brush). Daoud and Cotton predicted for star shaped polymers in the matrix a change in the scaling behavior as the blobs of the chains become non-ideal. The density profile is directly related to the average brush height, h, or the extension of the corona chains. Neglecting the contribution of the core to the radius of the star, they found that h ∝ N 1∕2 g 1∕4 b. Although the former relation exhibits "ideal" scaling with respect to the chain length dependence, the presence of the factor 1∕4 r shows that the radius is in fact larger than it would be for a single linear chain. Thus, although we are in a regime where the chain seems to be ideal, the structure is actually stretched. Vogiatzis and Theodorou [192] investigated the structural features of polystyrene brushes grafted on spherical silica nanoparticles immersed in polystyrene by means of a Monte Carlo methodology based on polymer mean field theory. The nanoparticle radii (either 8 or 13 nm) were held constant, while the grafting density and the lengths of grafted and matrix chains were varied systematically in a series of simulations. The primary objective of that work was to simulate realistic nanocomposite systems of specific chemistry at experimentally accessible length scales and study the structure and scaling of the grafted brush. In Fig. 12 the average thickness is plotted versus N 1∕2 g 1∕4 . N g is measured in Kuhn segments per chain and in chains per nm 2 . The scaling prediction of Daoud and Cotton seems to be fullfilled for both the rms height ⟨ h 2 ⟩ 1 2 and the height containing 99% of the brush material, h 99% . The dashed lines are linear fits, confirming the good agreement of the simulation data with the theoretical scaling behavior. The agreement seems to be better for the h 99% data points. This was expected, since the average brush thickness, as defined in Eq. (62), is more sensitive to the discretization of the model and to the post processing of the data, than the straightforward definition of the shell in which the 99% of the brush material can be found. Moreover the results for the h 99% estimate and the scattering pattern of the whole grafted corona were in favorable agreement with SANS measurements.
Voyiatzis et al. [193] studied the confinement induced effects on the accessible volume and the cavity size distribution in polystyrene-silica nanocomposites by atomistic Molecular Dynamics simulations. The composite systems contained a single -quartz silica nanoparticle, either bare or grafted with atactic PS chains, which was embedded into an unentangled atactic PS matrix [158]. Both free and grafted chains consisted of 20 monomers. The considered nanoparticle diameters were 3.0, 4.0 and 5.0 nm and three different grafting densities were studied: 0.0, 0.5 and 1.0 chains / nm 2 . Those authors investigated the cavity distribution size by employing four spherical probe particles. Apart from the limiting case of a dot-like probe particle (zero radius), the considered probe particles had radii, r p , of 0.128, 0.209 and 0.250 nm, corresponding to hard-sphere representations of helium, methane and ethane. The "unoccupied" volume was defined as the volume accessible to a probe particle of r p = 0.
The influence of the grafting density on the spatial distribution of the unoccupied volume fraction, v Un , and the specific volume, v Sp , for a nanoparticle with a diameter of 3 nm is presented in Fig. 13. The black horizontal line corresponds to the bulk value of v Un . The greatest reductions of v Un relative to the bulk value occur very close to the surface, at distances smaller than 1 nm. The variation Fig. 12 The calculated brush thickness (either ⟨ h 2 ⟩ 1 2 or h 99% ) is plotted versus the degree of polymerization of grafted chains, N g , times the grafting density, 1∕4 . Points correspond to systems containing an 8-nm-radius silica particle grafted with PS chains and dispersed in PS matrix. Linear behavior is predicted by the model proposed by Daoud and Cotton for star shaped polymers [181]. (Reprinted with permission from Ref. [192]. Copyright (2011) American Chemical Society.) in the v Un distribution of a grafted and an ungrafted nanoparticle is different. The separation from the particle for which v Un is below the bulk value for the grafted 1.0 chains/nm 2 system exceeds by approximately 50% the distance for the ungrafted system. The v Un profile for a grafting density of 0.5 chains / nm 2 lies in between the two extremes. Unlike v Un , the v Sp spatial distributions exhibit a strong dependence on the grafting density. An increase of the grafting density leads to increased v Un values close to the particle's surface. This behavior was attributed to the (i) the chemistry of the employed linker molecule and (ii) the expulsion of the free chains from adsorbing on the nanoparticle surface. Contrary to intuitive expectations, variations of the accessible volume were not directly related to changes of the specific volume.

Dynamics
A complete understanding of PNC dynamics requires confronting the difficult many-body problem associated with non-dilute particle concentrations and coupled nanoparticle and polymer motions over many time-and length-scales [194][195][196]. Simulations are a valuable option, but are computationally very intensive, resulting in only a limited parameter range being tractable to study, typically involving rather small particles and weakly entangled polymers. The transport properties of nanoparticle-polymer mixtures have been the focus of much recent attention [197][198][199][200][201][202][203][204][205][206][207]. Central problems in the area are the diffusion of nanoparticles and polymers through the nanocomposite melts, as well as the local polymer dynamics in the proximity of the filler particles. For example, the incorporation of nanoparticles can strongly modify the viscosity of polymer melts [208], and the center-of-mass mobility of polymer chains can be strongly retarded depending on nanoparticle size and concentration [209,210].

Experimental Findings
There is good understanding of the motion of very large or very small colloidal particles of radius R n in a polymer melt. The nanoparticle diffusion coefficient, D, in the large particle limit follows the classic continuum Stokes-Einstein relation [211]. For a large and massive solute molecule of radius R n in a solvent consisting of much smaller and lighter molecules, the diffusion coefficient, D, of the solute is given by [212]: where k B is Boltzmann's constant, T is the absolute temperature, is the solvent viscosity and f takes the values of 4 or 6 for slip or stick boundary conditions at the solute surface, respectively [203]. The corresponding behavior of small nanoparticles, comparable to the size of a monomer, is also described by the Stokes-Einstein relationship but with a length-scale dependent viscosity that is smaller than the macroscopic bulk value [213]. The relevant apparent viscosity is controlled by the relaxation of subsections of chains with an end-to-end distance comparable to the nanoparticle size, as has been verified by Molecular Dynamics simulations [214]. Understanding nanoparticle diffusion in the polymer matrix is of fundamental importance.
Despite the rather good understanding of the diffusion of the particles in the two limits (very large and very small), the dynamical behavior of nanoparticles of size comparable to the entanglement mesh size of the polymer is contentious [205,[215][216][217]. Brochard-Wyart and de Gennes [213] argued that the particle diffusion constant follows normal Stokes-Einstein behavior essentially when its size becomes larger than the entanglement mesh size. Such a sharp sizedependent crossover to Stokes-Einstein scaling has been observed by Szymanski et al. [218]. On the contrary, Cai et al. [201] speculated that the motion of these intermediate sized nanoparticles should be faster than Stokes-Einstein behavior, since diffusion can be facilitated by hoplike motions trough the polymer's entanglement mesh. The latter is also supported by a microscopic force-level theory, 13 Distribution of unoccupied volume fraction, v Un , (lefty-axis) and the specific volume, v Sp , (righty-axis) in the vicinity of a nanoparticle with a diameter of 3 nm at a termperature of 590 K. Thehorizontal line corresponds to the PS bulk value of v Un . The grafting density is varying from 0.0 (ungrafted) to 0.5 and 1.0 grafted chains / nm 2 . (Reprinted from [193] with permission from Elsevier.) wherein chain relaxation and local entanglement mesh fluctuations, i.e. "constraint release", dominate over hopping [219].
Somoza et al. [220] studied experimentally the anthracene rotation in poly(dimethylosiloxane) and poly(isobutylene) by gradually increasing the chain length. These authors reported that the diffusivity of the particles exhibited a sharp transition with the increase of the polymer radius of gyration, R g , becoming dependent on the "nanoviscosity" (rotation time of dissolved anthracene was used as a measure of the viscosity on a nanometer-sized object) rather than the macroscopic viscosity for small R n ∕R g ratios (with R n being the particle size). Narayaman et al. [221] used X-ray photon correlation spectroscopy in conjunction with resonance-enhanced grazing-incidence small-angle X-ray scattering to probe the particle dynamics in thin films, and also found that the particle dynamics differ from the Stokes-Einstein Brownian motion, the difference being caused by the viscoelastic effects and interparticle interactions. Meanwhile, Tuteja et al. [198] reported that the nanoparticles diffuse two orders of magnitude faster in a polymer liquid than the prediction of the Stokes-Einstein relation, an observation possibly attributable to the nanoparticles being smaller than the entanglement mesh. Later, Grabowski et al. [199] also observed strong enhancement (250 times) of the diffusion of gold nanoparticles in poly(butyl methacrylate), under conditions where the nanoparticle dimensions were smaller than the entanglement mesh length of the polymer.
Cai et al. [201] have carried out an extensive study of nanoparticle diffusion by employing scaling theory to predict the motion of a probe nanoparticle of size R n experiencing thermal motion in polymer solutions and melts. Particles with size smaller than the solution correlation length, , undergo ordinary diffusion with a diffusion coefficient similar to that in pure solvent. The motion of particles of intermediate size ( < d < pp ), where pp is the tube diameter for entangled polymer liquids, is subdiffusive at short time scales, since their motion is affected by subsections of polymer chains. At long time scales the motion of these particles is diffusive, and their diffusion coefficient is determined by the effective viscosity of a polymer liquid with chains of size comparable to the particle diameter R n . The motion of particles larger than the tube diameter pp at time scales shorter than the relaxation time e of an entanglement strand is similar to the motion of particles of intermediate size. At longer time scales (t > e ) large particles (d > pp ) are trapped by the entanglement mesh, and to move further they have to wait for the surrounding polymer chains to relax at the reptation time scale rep . At longer times t > rep , the motion of such large particles (d > pp ) is diffusive with diffusion coefficient determined by the bulk viscosity of the entangled polymer liquids. Finally, for nanoparticles with diameters larger than the entanglement mesh size it appears that the competition of full chain relaxation versus the nanoparticle hopping through entanglement gates controls nanoparticle diffusion [222].

Insight Obtained from Simulations
Liu et al. [214] employed Molecular Dynamics simulations of the Kremer-Grest model [149] in order to investigate the diffusion process of spherical nanoparticles in polymer melts. Their results indicated that the radius of gyration of the polymer chains was the key factor determining the validity of the Stokes-Einstein relation in describing the particle diffusion at infinite dilution. In Fig. 14 the diffusion coefficient estimated by the MD simulations is presented alongside the predictions of the Stokes-Einstein formula. It was found that, with the increase of the size ratio of R n ∕R g , the Stokes-Einstein diffusion coefficient gradually approximates the MD data under the slip (dotted curve) boundary condition. The use of purely repulsive non-bonded interactions fully justifies the use of the slip ( f = 4), instead of the stick ( f = 6), boundary condition. As the size ratio, R n ∕R g increases to 1, the predicted diffusion coefficients agree reasonably well with those extracted from the simulations. However, in the small size ratio, large deviation is observed which qualitatively agrees with the experimental results [198]. It seems like the local viscosity experienced by nanoparticles is much smaller than the macroscopic viscosity, as speculated by Wyart and de Gennes [213] and other researchers [223][224][225][226]. Finally, it should be noted that in the experiments chains are strongly entangled, in Fig. 14 The diffusion coefficient, D, of a spherical particle as a function of the ratio R n ∕R g . The particle mass is proportional to its volume. Open squares represent Molecular Dynamics data, while full dots represent the Stokes-Einstein relation predictions with slip boundary condition. (Reprinted with permission from Ref. [214]. Copyright (2008) American Chemical Society.) the reptation regime, while the polymer length used by Liu et al. is smaller than the entanglement length of the polymer chain [227].
Yamamoto and Schweizer [219,228] have formulated and applied a microscopic statistical-mechanical theory, based on the Polymer Reference Interaction Site Model (PRISM) integral equation theory [229], for the nonhydrodynamic relative diffusion coefficient of a pair of spherical nanoparticles in entangled polymer melts. Their work was based on a combination of Brownian motion, mode-coupling, and polymer physics ideas. They focused on the mesoscopic regime, where particles are larger than the entanglement spacing. The overall magnitude of the relative diffusivity was controlled by the ratio of the particle to tube diameter and the number of entanglements per chain. Figure 15 presents model calculations of the total relative diffusivity, D (rel) (h) as a function of h∕2R n (with h being the interparticle surface-to-surface separation distance) for two reduced particle diameters. The ordinate of the figure is normalized by the single particle Stokes-Einstein result, D SE , while the abscissa extends up to the point where the theory is argued to be reasonable. That theory is based on the mode-coupling idea that the relevant slow dynamical variable is the bilinear coupling of the nanoparticle and the collective polymer density fluctuations. The original approach [219] was not self-consistent since it assumed that the constraining forces on a particle relax entirely due to the lengthscale dependent motions of the polymer melt (constraint release regime), which is an accurate simplification when particles are larger than d T . Figure 15 exhibits several interesting trends. First, the relative diffusivity approaches the asymptotic value D (self) ∕D SE at h∕2R n >> 1, verifying the "isolated" limit (two particles at infinite dilution). Note that it does not necessarily approach unity if a Stokes-Einstein violation is present at the single-particle level. The overall deviation of D (rel) from the hydrodynamic result is enhanced as 2R n ∕d T decreases or N∕N e increases, in analogy with a singleparticle Stokes-Einstein violation effect. The underlying physical mechanism can be understood as small nanoparticles acquiring high mobility due to a weaker coupling to the slow relaxation of entangled melts compared to the continuum theory. As a consequence, the overestimate of friction by a hydrodynamic approach grows as particle size decreases and/or chain length increases. The second general feature in Fig. 15 concerns the role of the number of entanglements, N∕N e . Deviations cannot be discerned from the hydrodynamic behavior for weakly entangled cases (N∕N e ∼ 1) for either particle size. One may physically rationalize this result by recalling that the Rouse-like collective relaxation is diffusive above the segmental lengthscale. Finally, the most important feature of Fig. 15 is the predicted non-trivial mobility enhancement compared to the hydrodynamic result over a wide range of h∕2R n . Hydrodynamics predicts zero relative diffusivity as h → 0 , whereas the results of Yamamoto and Schweizer remain non-zero down to small separations.
Kalathi et al. [230] have employed large-scale molecular dynamics simulations in order to examine the role of entanglements on nanoparticle dynamics in the crossover regime, where the diameter of the particles, NP , is larger than 2 − 10d T with d T being the entanglement tube diameter. The transport behavior of nanoparticles in the crossover size limit appears to be complicated by hopping effects, length-scale dependent entanglement forces and dynamics, and the interactions of polymers and nanoparticles. These authors simulated weakly interacting mixtures of nanoparticles and bead-spring polymer melts. For the polymer melts considered in that work, the entanglement chain length is approximately 45, N e~4 5, and d T in units of monomer diameter, , is around 7 (the nanoparticle diameters were NP = 1 − 15, in units of polymer ). The diffusion coefficients of nanoparticles smaller than d T ∼7 − 10 (i.e. NP = 1, 3, and 5, respectively, Fig. 16(a)) in long chain melts show that the relevant viscosity corresponds to a section of the chain with N NP monomers that satisfies 2 NP = N NP 2 . For shorter chains the data can be described where = 1 N (where 1 is the viscosity of a monomer fluid at the same density) and f ∼4 (asymptote in Fig. 16(a)). The results of Kalathi et al. [230] for smaller nanoparticles are in good agreement with the predictions of the theory of Yamamoto and Schweizer [228]. Figure 16(b) presents the results for nanoparticles with sizes larger than the entanglement mesh length, d T ∼7 − 10. The diffusion of these particles in the longer chain melts does not follow the "universal" plateau seen for small nanoparticles. These results suggest that the chain-scale dynamics does not control nanoparticle diffusion (no Stokes-Einstein scaling). However, the fact that the diffusivity at high N of these intermediate-sized particles does not reach the same plateau as the small nanoparticles (Fig. 16) suggests that another effect, probably entanglements, plays an important role. Despite the fact that no conclusive evidence of hopping-controlled transport was found, the spontaneous fluctuations of the entanglement mesh (constraint release) in the moderately long chain melts, may be the most important mode of nanoparticle transport through the bulk of the material, in agreement with theoretical predictions [228].

Polymer Diffusion and Dynamics
Early theoretical studies [231][232][233] have shown that polymer diffusion through heterogeneous media is slowed down due to entropy losses associated with impenetrable obstacles. A natural choice of a parameter for quantifying the ability of a nanocomposite to slow down polymer diffusion is the spacing between the surfaces of neighboring nanoparticles [234]. For monodisperse nanoparticles, this spacing can be defined as the interparticle distance, d inter , given by [235]: where d and are the nanoparticle diameter and nanoparticle volume fraction, respectively.The maximum packing density of the nanoparticles, max , depends on the packing type, such as simple cubic ( max = 0.524), face-centered cubic ( max = 0.740), body-centered cubic ( max = 0.680), and random dense packing ( max = 0.637). Independently of max , d inter decreases as nanoparticle size decreases at fixed , suggesting that smaller nanoparticles slow down polymer diffusion more effectively than larger ones [231,232].

Experimental Findings
Gam et al. [210] have measured the tracer diffusion of deuterated polystyrene (dPS) in a polystyrene nanocomposite containing silica nanoparticles, with number average diameters, d n , of 28.8 and 12.8 nm, using elastic recoil detection. The corresponding volume fractions of the large and small nanoparticles, , ranged from 0 to 0.5, and 0 to 0.1, respectively. At the same volume fraction of nanoparticles, the tracer diffusion of dPS is reduced as nanoparticle size NP should be a constant, independent of chain length [201,214]. a Nanoparticles that are smaller than the entanglement mesh size, Nanoparticles that are larger than the entanglement mesh size, NP ≥ d T . In all cases, a slip boundary condition was assumed. Predictions of the theory of Yamamoto and Schweizer [228] are presented in solid lines. Moreover, if we estimate the Rouse viscosity based on the actual chain length N, we get D * decreases because the interparticle distance between nanoparticles, d inter , decreases. The reduced diffusion coefficient, defined as the tracer diffusion coefficient in the nanocomposite relative to pure PS (D∕D 0 ) is plotted against the confinement parameter divided by the tracer size in Fig. 17. All measurements nearly collapse onto a master curve [234], although D∕D 0 is slightly higher for the smaller particles. For d inter = ID < 2R g , D∕D 0 decreases rapidly as ID∕2R g decreases. For ID > 2R g , D∕D 0 remains less than 1 indicating that entropy loss reduces diffusion even when nanoparticles are far apart relative to the tracer size. The dashed line is an empirical fit, because a theory relating D to the fundamental system parameters is lacking.
Schneider et al. [236] studied experimentally the relaxation of entangled poly(ethylene-alt-propylene) (PEP) chains (tube diameter ~5 nm) filled with silica nanoparticles (average diameter ~17 nm). The silica volume fraction was varied between 0.0 and 0.6 (as that was estimated from the measured weight fraction of silica in the nanocomposite). Neutron spin echo spectroscopy (NSE) was empoloyed in order to explore chain dynamics in these nanocomposites, characterized by non-attractive interactions. The resulting collective dynamic scattering function data were analyzed by employing the idea of a tube-like confinement for chain relaxation below the reptation time. The following conclusions were drawn from their study: (i) the monomeric relaxation rates were not unaffected by the addition of nanoparticles, even at high particle loadings; (ii) chain conformations remain Gaussian for all loadings considered; and (iii) the tube diameter determined from analysis of neutron spin echo data decreases monotonically upon adding nanoparticles. Two contributions to overall chain dynamics were speculated. On the one hand, the number of topological chain-chain entanglements decreases with increased nanoparticle loading, i.e., the chains disentangle from each other since a part of the system volume is occupied by the NPs. On the other hand, the chain acceleration caused by the reduction of entanglements is (more than) compensated by the geometric constraints that nanoparticles present to chain dynamics. Since that second factor dominates at large loadings, the neutron scattering experiments suggested an increase in chain relaxation time, while at the same time a reduction of chain-chain entanglements and an increase of particle-chain entanglements take place.

Insight Obtained from Simulations
Desai et al. [237] investigated the chain dynamics of Kremer-Grest polymer melts, composed of chains with a relatively high degree of polymerization (N = 80) filled with solid nanoparticles using molecular dynamics simulations. These authors found that chain diffusivity is enhanced relative to its bulk value when polymer-particle interactions are repulsive and is reduced when polymer-particle interactions are strongly attractive (Fig. 18). In both cases chain diffusivity assumes its bulk value when the chain center of mass is about one radius of gyration, R g , away from the particle surface. As shown in Fig. 18 for the particle volume fraction of 10%, the average chain diffusion coefficient is reduced by a factor of 2 in the presence of strongly attractive particles. The case of repulsive particles appears to be even more interesting, where the diffusion coefficient initially increases with increasing particle concentration, but then reaches a maximum before decreasing with further increase of the particle concentration. While the initial particle concentration dependence of the diffusion coefficient reflects the polymer-particle interactions, higher concentrations always lead to a reduction of the diffusion coefficient, which may be attributed to geometrical reasons (i.e. the presence of tortuous paths in systems with high particle loadings).
Kalathi et al. [239] employed large-scale molecular dynamics simulations in order to study the internal relaxations of chains in nanoparticle/polymer composites. They examined the Rouse modes of the chains, which resemble the observables of the self-intermediate scattering function, typically determined in an (incoherent) inelastic neutron scattering experiment. The Rouse modes, p = 0, 1, 2, ..., N − 1, of a chain of length N are defined as [240]: The time autocorrelation of the Rouse modes is predicted to decay exponentially and independently for each node p for an ideal chain, with relaxation time, p . The p = 0 mode describes the motion of the chain center-of-mass, while the modes with p ≥ 1 describe internal relaxations with a mode number p corresponding to a sub-chain of (N − 1)∕p segments. The comparison of the relaxation times of the different modes for chains in the PNCs for three different degrees of polymerization, N, filled with nanoparticles of different sizes for = 0.1 to neat melt is presented in Fig. 19(a)-(c). Their results (Fig. 19) showed that, for weakly interacting mixtures of nanoparticles and polymers, the effective monomeric relaxation rates are faster than in neat melt when the nanoparticles are smaller that the entanglement mesh size. In this case, the nanoparticles serve to reduce both the monomeric friction and the entanglements in the polymer melt, as in the case of polymer-solvent mixtures. On the contrary, for nanoparticles larger than half the entanglement mesh size, the effective monomer relaxation remains unaffected for low nanoparticle concentrations. Even in this case, strong reduction of chain entanglements was observed. These authors concluded that the role of nanoparticles is to always reduce the number of entanglements. By assuming that the relaxation time for a chain follows the crossover bridging Rouse to reptation dynamics, the large p modes directly yield information on the monomer friction and in the limit of p = 1, the plateau of eff p ∕ eff p,neat is directly proportional to the ratio of 0 ∕N e in the PNC compared to that in the pure melt. For small nanoparticles, which act as a diluent, there was an additional speedup, which was attributed to a reduction in entanglements, quantified by N e, melt ∕N e, PNC ∼0.9 for long chains, which can also be extracted by the stretching exponents (Fig. 19(e)).

Insight Obtained from Simulations
Brown et al. [167,241] were among the first to study the local dynamics of a model nanocomposite system. They examined the structure and dynamics of a system containing an inorganic (silica) nanoparticle embedded in a polymer (polyethylene-like) matrix. They thoroughly discussed the variation of structure and dynamics with increasing distance from the polymer-particle interface and as a function of pressure. A clear structuring of the linear polymer chains around the silica nanoparticle was observed, with prominent first and second peaks in the radial density function and concurrent development of preferred chain orientation. Evidence of chain immobilization was less obvious overall, although dynamic properties were more sensitive to changes in the pressure. Long simulations were carried out to determine the variation in the glass transition of the filled polymers as compared to the pure systems. In Fig. 20 the average relaxation times of the torsional autocorrelation function are presented, resolved into concentric shells of 5 Å thickness around the nanoparticle center of mass. This assignment was based on the position of the center of mass of the four atoms involved in the torsion at the time origin of the observation. Although this means that at some later time an angle may belong to a different shell, it avoids the bias that would result from selecting only angles that remain in a particular shell (in any case diffusion of chains is relatively slow so that should not be a problem). Based on Fig. 20, there is some indication that the decreased translational mobility near the interface increases the relaxation time associated with torsional equilibration (establishment of the trans-gauche equilibrium in the systems containing the nanoparticles with diameters of 3 and 6 nm (R30L and R60L, respectively), but otherwise most of the characteristic times are very close to those obtained for the neat system, in agreement with previous studies [241,242] of the same authors. It was concluded that, within errors, the interphase thickness was independent of the size of the nanoparticle for the range of particle systems analyzed.
Vogiatzis and Theodorou [156] produced atomistic configurations of fullerene-filled polystyrene melts by reverse mapping well-equilibrated coarse-grained melt configurations, sampled by connectivity altering Monte Carlo, to the atomistic level via a rigorous quasi Metropolis reconstruction. The main goal of their work, i.e., the study of PS-C 60 dynamics at the segmental and local levels, has been accomplished by analyzing long MD trajectories of their well-equilibrated reverse-mapped structures. Their simulation results generally indicate that the addition of C 60 to PS leads to slower segmental dynamics (as estimated by characteristic times extracted from the decay of orientational time-autocorrelation functions of suitably chosen vectors), suggesting an increase of the bulk glass transition temperature, T g , by about 1 K, upon the addition of C 60 at a concentration of 1% by weight (in favorable agreement with differential scanning calorimetry measurements [243]). They then employed a space discretization in order to study the effect of C 60 on segmental dynamics at a local level. Each fullerene served as the center of a Voronoi cell, whose volume was determined by the neighboring fullerenes. Backbone carbons lying in every cell were tracked throughout the atomistic Molecular Dynamics trajectory and their mean-square displacement (MSD) was measured for the time they resided in the same cell. Overall mean-square displacement of backbone atoms was found to be smaller in the presence of fullerenes, than in bulk PS. However, atoms moving in smaller (more confined) Voronoi cells exhibited faster motion than the atoms moving inside larger Voronoi cells. Figure 21 presents the MSD of backbone carbon atoms as a function of time at a temperature of 480 K for both the filled and unfilled systems. As can be seen, nanocomposite systems exhibit lower mobility when compared to their neat counterparts. The MSD of backbone carbons is depressed upon the addition of fullerenes, in good agreement with the neutron scattering observations of Kropka et al. [243]. In the inset to Fig. 21, a logarithmic plot of the MSD is presented. The scaling of t 1∕2 is expected for the very short time behavior studied, as Likhtman and McLeish [244] have estimated that the time marking the onset of the effect of topological constraints on segmental motion, e , is 3.36 × 10 4 s for polystyrene.
Moreover, Vogiatzis and Theodorou [156] estimated the local mean-square displacement (MSD) of backbone carbon atoms of PS, for the timespan an atom spends inside a particular cell of the Voronoi tessellation. In their analysis they used the average MSD from the three most confined and three least confined cells. They observed that the volume of the Voronoi cells did not change significantly as a function of time. Based on that analysis for the nanocomposite system, the degree of depression was found to be a function of the confinement induced by the fullerenes. The diffusion of chains was spatially inhomogeneous, as observed by Desai et al. [237] earlier. Small Voronoi cells lead to higher mobility of the polymer segments within them. Despite the fact that the addition of fullerenes limited the diffusion of polymeric chains, there existed regions in space, where the polymer could recover part of its dynamics due to the high level of confinement. This finding was then correlated with the increased rotational diffusion of fullerenes, as the volume of the Voronoi cells became smaller. These authors envisioned fullerenes as nanoscopic millstones, forcing the polymeric chains to diffuse faster, when close to them, while the geometrical constraints imposed by their presence force the chains to diffuse more slowly.
Pandey et al. [245] have extensively studied the local dynamics and the conformational properties of polyisoprene next to a smooth graphite surface constructed by graphene layers, via a multiscale simulation methodology. These authors first performed fully atomistic molecular dynamics simulations of isoprene oligomers, next to the surface. Subsequently, Monte Carlo simulations of a systematically derived coarse-grained model were employed in order to create several uncorrelated structures for polymer systems. A reverse mapping strategy was developed in order to reintroduce atomistic detail into the coarsegrained configurations. Finally, multiple extensive fully atomistic simulations with large systems of long macromolecules were conducted to examine local dynamics in Fig. 20 The radial dependence of the relaxation time of the torsional autocorrelation function of polyethylene around a silica nanoparticle. All values are averages taken in 5 Åshells around the nanoparticle center of mass. The points have been offset slightly for the three systems along the x axis for clarity. The dotted line simply indicates the value obtained for the 30-chain pure polymer system at the same temperature (400 K). (Reprinted with permission from Ref. [167]. Copyright (2008) American Chemical Society.)

Fig. 21
Mean-squared atomic displacements (MSD) of backbone carbon atoms as a function of time for filled and unfilled polystyrene systems at T = 480 K. In the case of fullerene nanocomposites, an analysis of the dependence of backbone MSD on confinement is also presented for most and least confined Voronoi cells (indicative error bars also included). In the inset to the figure, the same data are presented in logarithmic axes. (Reprinted with permission from Ref. [156]. Copyright (2014) American Chemical Society.) proximity to graphite. Their findings supported the presence of increased dynamic heterogeneity emerging from both intermolecular interactions with the flat surface and intramolecular cooperativity. For each system, Pandey et al. extracted bond orientation autocorrelation functions and sorted them in intervals of 0.03 nm based on the position of the midpoint of the c-CH bond throughout the simulation trajectory. For each interval, the autocorrelation curves were averaged weighting by the population of c-CH vectors found in a specific interval from each run. The mean correlation times increased substantially in the proximity of the surface, with dynamics at the surface almost 20 times slower (independently of the molecular weight of the chains) than in bulk PI. Figure 22 presents a qualitative visual inspection of the distribution of times for a specific configuration. Pandey et al. [245] evaluated an autocorrelation function individually for every c-CH vector and colored each segment from blue to red signifying lower relaxation times and higher mobility. As shown, segments in proximity to the surface were found to be slower. They then separated all train segments based on their length and calculated correlation times for each position along the length of the train segment (Fig 22(b)), which was an extremely challenging task. Despite a significant statistical error, several features are evident in Fig.22(b). Specifically, when a chain makes a single contact, dynamics is only decelerated to a small extent. The second important finding was that, as train segments grow in length along the surface, the dynamics of the repeat units becomes progressively slower, with the findings implying that, similar to chain-ends in bulk dynamics, [246,247] the ends of train segments contribute to increased dynamic heterogeneities on the surface. However, the former are only significant for short chains, the latter were present for any chain length studied. In addition, short train segments can be more pronounced around surfaces with higher curvature [16,159]. Finally, a PI specific result was the asymmetry present along a train segment originating from the methyl group, much alike findings on bulk dynamics along the chain backbone [247].
Rissanou and Harmandaris [248] presented a detailed analysis of the dynamics of three different polymer-graphene systems, through atomistic Molecular Dynamics simulations. In more detail they studied (a) PS-graphene, (b) PMMA-graphene and (c) PE-graphene interfacial systems, as well as the corresponding bulk polymer systems. For PS and PMMA polymer chains were 10-mers while PE chains were 20-mer, in order to ensure that the backbone consisted of almost the same number of CH 2 , and/or CH groups for all systems (i.e. approximately 20 in all cases). A characteristic quantity of the molecular level is the endto-end vector ee (t), whose autocorrelation function provides information for the orientational dynamics at the entire chain level. Rissanou and Harmandaris performed an analysis of end-to-end autocorrelation function at different adsorption layers and fit the corresponding curves for all chains to the Kohlrausch-Williams-Watts (KWW) function [249][250][251]. At the entire chain level, the integral below the KWW curves defines the molecular chain endto-end relaxation time, ee mol . The molecular relaxation times together with the stretching exponent, , of the KWW fits are presented in Fig. 23 (a) and (b) as functions of the distance from the surface. Data in Fig. 23(a) reveal the dramatic increase of ee mol close to the graphene layer, compared to the corresponding bulk values, shown with dashed lines. Furthermore, a slight difference in the distance at which the ee mol reaches the plateau distance-independent bulk value was observed: for PE it is about ∼2 nm, whereas for PMMA and PS is about 3 − 4 nm. The extreme difference in relaxation times between PE and the other two systems is obvious. The exponent of the KWW relation for PE and PS reaches a constant value in the bulk region, while the same does not apply for PMMA. These values show that in the bulk region PMMA has the wider distribution of relaxation times, PS follows and PE has the narrowest one.

Phase Behavior
Polymer-nanoparticle blends exhibit a rich phase behavior which is directly tied to the thermal, mechanical, and optical properties of the composite system, with the achievement of uniform dispersion being a long-standing challenge [1,14,15,33,150,252,253].
Significant progress towards the development of microscopic predictive theories of the equilibrium structure and phase behavior of polymer nanocomposites has been made recently based on liquid state integral equation formulations, density functional calculations and self-consistent mean field approaches. All these methods can complement or surpass the explicit atom methods like Monte Carlo and Molecular Dynamics, which have the potential to quantitatively predict structural correlations, thermodynamics and phase behavior.
Chatterjee and Schweizer were the first to develop an analytical integral equation theory for treating polymerinduced effects on the structure and thermodynamics of dilute suspensions of hard spheres [254]. Results were presented for the potential of mean force, free energy of insertion per particle into a polymer solution, and the second virial coefficient between spheres. Later, Hoopper et al. [255] employed the Polymer Reference Interaction Site Model (PRISM) theory to investigate structure, effective forces, and thermodynamics in dense polymer-particle mixtures in the one and two particle limit [144,145].

Bare Nanoparticles
Hall et al. [256,257] employed Polymer Reference Interaction Site Model (PRISM) liquid state theory to study phase transitions and structure of dense mixtures of hard nanoparticles and flexible polymer coils. Their calculations were performed over the entire compositional range from the polymer melt to the hard sphere fluid, with the focus being on polymers that adsorb on nanoparticles. Many body correlation effects were fully accounted for in the determination of the spinodal phase separation instabilities. An example phase diagram is presented in Fig. 24. It can be discerned that depletion and bridging phase separation occur at low and high attraction strengths, respectively. Quantitatively, many particle effects are found to always reduce miscibility. Depletion phase separation was similar for different attraction ranges, with a critical point at rather low filler volume fractions. This is in contrast to the bridging induced demixing transition where the critical point is located at very high nanoparticle volume fractions. Moreover, increasing the attraction range increases the thickness of the bound layer and the importance of many body effects, which further decreases miscibility in the high filler volume fraction regime relative to what was predicted by a two particle virial analysis [145]. However, when bridging effects are very strong and phase separation occurs at low volume fractions, decreasing the attraction range can lead to a stronger, shorter range bridging attraction that reduces miscibility. Increasing particle size generally disfavors miscibility on both the depletion and bridging sides of the spinodal phase diagram, though the effect on depletion is more significant.
Wei et al. [258] have investigated silica nanoparticle dispersions in polystyrene, poly(methyl methacrylate), and poly(ethylene oxide) melts by means of a density functional approach. The polymer chains were regarded as coarse-grained semi-flexible coils whose segment size matched the Kuhn length of the polymer under investigation. The particle-particle and particle-polymer interactions were calculated in the grounds of the Hamaker theory, following Vogiatzis and Theodorou [148,192]. In order to characterize nanoparticle dispersion, Wei et al. employed the second virial coefficient, B 2 , defined as: where the first term accounts for the particle contribution, and the second one is the polymer mediated contribution. The local density of particles is denoted as n ( ), while the average particle density as n . n ( ) varies as a (67) function of distance between particles, making it a critical link to particle microstructure. The pair correlation function approaches asymptotically n ( )∕ n = 0 when r < 2R n (R n being the radius of the particles) as particles cannot interpenetrate and n ( )∕ n ≃ 1 as r → ∞ as the likelihood of finding a particle becomes proportional to the average particle density. Figure 25 shows the second virial coefficients for different particle sizes at constant particle volume fraction = 5%. Positive values of B 2 indicate stable particle dispersion (effective particle-particle repulsion), while a negative value signifies unstable dispersion. The results are in agreement with previous theoretical studies: the tendency to dispersion increases as the particle size increases [259]. It can be seen that B 2 becomes independent of the particle size after a threshold value, that meaning the effect of particle size on the pair correlation function becomes insignificant. Before that critical value, B 2 increases but its increasing amplitude declines as the particle size increases. Density functional theory confirms that large particles are more likely to achieve stable dispersion than small particles.

Polymer Grafted Nanoparticles
One approach for controlling the particle dispersion in the polymer matrix is to alter the particle surface chemistry through the attachment of polymer chains. The composition, architecture, and distribution of the grafted chains can be carefully designed to tailor interparticle interactions, thereby controlling the dispersion state [33]. In the special case where the chemical composition of the graft and matrix chains are identical, the entropic contributions dominate the thermodynamics [175,260], and uniform dispersion can be achieved both for the case of spherical nanoparticles [174,261], and the case of nanorods [262].

Experimental Findings
In a number of studies, empirical phase diagrams have been developed for this special case, where the particle miscibility is a function of grafting density, , the ratio between the molecular weight of the grafted chains, N, and matrix chains, P, i.e. P / N, as well as the particle radius, R n [14,71,72,173,174,261,[263][264][265][266]. Sunday et al. [174,267] quantified the stability of polystyrene-grafted silica nanoparticles in PS matrices with ultrasmall angle X-ray scattering (USAXS) and transmission electron microscopy (TEM). They developed the phase diagram presented in Fig. 26 to predict nanoparticle dispersion based on the graft polymer density, , and the graft and free polymer molecular weights, or N and P, respectively. The phase diagram of Fig. 26 shows three distinct regions. When is below the allophobic limit ( 1 ), the surface coverage of grafted chains is low enough that the interactions between the particle and the matrix chains are not screened out sufficiently and the grafted particles behave similarly to block copolymers, aggregating into distinct morphologies [14,268,269]. As increases above 1 , the matrix and grafted chains interpenetrate, resulting in a "wet" brush and repulsive interactions between nanoparticles, thus stabilizing their dispersion. The autophobic dewetting line corresponds to a continuous, second-order transition, resulting from the expulsion of the melt from the brush for densely grafted chains ( N 1∕2 > 1) which should lead to nanoparticle aggregation through the attraction between graft layers [270]. Larger values of R n , , or P/N result in a larger entropic penalty for intermixing due to crowding of the grafted layer. As the entropic penalty grows, the interpenetration width between the matrix and grafted chains decreases until the matrix chains are completely expelled, resulting in attractive interactions and particle aggregation. This occurs above the autophobic phase transition at 2 , a discontinuous, first-order transition at low grafting densities.
Bansal et al. [27] have experimentally observed and modeled the anisotropic self-assembly of small PS-gsilica nanoparticles with cores of R n ∼10 − 13nm in the allophobic dewetting region at lower grafting densities ( = 0.01 − 0.10 chains∕nm 2 ) [14]. Using slightly higher grafting densities ( ∼0.2 − 0.7 chains∕nm 2 ), Chevigny et al. [186] used similar-sized PS-g-silica NPs with N = 5 − 50 kg/mol in P = 140 kg/mol where particles with the longest grafts (P∕N = 2.8) dispersed uniformly, whereas those with the shortest grafts (P∕N = 28) phase separated from the bulk, forming spherical aggregates. The dispersion of silica NPs with higher graft densities has been investigated in which two sets of PS-g-silica NPs, the first with N = 110 kg/mol and = 0.27 chains∕nm 2 and the second with N = 160 kg/mol and = 0.57 chains∕nm 2 have been shown to disperse at least up to P∕N = 2.3 [271] and 1.6 [27], respectively.
Another way of tuning the mechanical properties of composite materials is by dispersing hydrophilic nanofillers in highly hydrophobic polymer matrices [272]. Martin et al. [273] have performed simulations and experiments on mixtures containing polymer grafted nanoparticles in a chemically distinct polymer matrix, where the graft and matrix polymers exhibit attractive enthalpic interactions at low temperatures that become progressively repulsive as temperature is increased.

Insight Obtained from Simulations
Trombly and Ganesan [264] have calculated the potential of mean force (PMF) between grafted nanoparticles immersed in a chemically identical polymer melt using a numerical implementation of polymer mean-field theory. These authors focused on the interpenetration width between Fig. 26 Illustration of the phase diagram for nanoparticle stability as a function of grafting density ( ) and the ratio of the lengths of the free over the grafted chains (P/N). Particles at low grafting densities encounter the allophobic dewetting transition at 1 . Increasing leads to complete wetting of the brush by the melt, stabilizing the nanoparticle dispersion. Increasing the graft density further leads to the autophobic dewetting transition at 2 and unstable dispersion. Particle dispersion is unstable at all grafting densities when P∕N > (P∕N) * . (Reprinted with permission from Ref. [174]. Copyright (2012) American Chemical Society.) the grafted and free chains and its relationship to the polymer-mediated interparticle interactions. To this end, they quantified the interpenetration width as a function of particle curvature, grafting density, and the relative molecular weights of the grafted and free chains.
Meng et al. [266] used Molecular Dynamics simulations to delineate the separation dependent forces between two polymer-grafted nanoparticles in a homopolymer melt and the associated potential of mean force (PMF). The nanoparticle radius (=5 in units of the chain monomers) and grafted brush length (=10) were held constant, while the grafting density and the polymer matrix length were varied systematically in a series of simulations. At first, it was shown that simulations of a single nanoparticle did not reveal any signatures of the expected autophobic dewetting of the brush with increasing polymer matrix length (in agreement with Monte Carlo simulations of Vogiatzis and Theodorou [192]). In fact, density distributions of the matrix and grafted chains around a single nanoparticle appeared to only depend on the grafting density, but not on the matrix chain length. However, the calculated forces between two nanoparticles in a melt, presented in Fig. 27, showed that increasing the matrix chain length, M, from 10 to 70 causes the interparticle PMF to go from purely repulsive to attractive with a well depth of the order of k B T (with k B being the Boltzmann constant). It was speculated that these results were purely entropic in origin and arise from a competition between brush-brush repulsion and an attractive inter-particle interaction caused by matrix depletion from the inter-nanoparticle zone (i.e. an Asakura-Oosawa type inter-particle attraction). Figure 27 compares the PMF with the results from Smith and Bedrov' simulations [274] of a similar coarse-grained system using the umbrella sampling method for an apparently identical chain length and coverage. The two studies are in qualitative agreement to each other, with the PMF of Meng et al. consistently shifted toward smaller separations.
Martin et al. [275] presented an integrated theory and simulation study of polydisperse polymer grafted nanoparticles in a polymer matrix to demonstrate the effect of polydisperisty in graft length on the potential of mean force between the grafted particles. It is evident from Fig. 28 that increasing polydispersity in graft length reduces the strength of repulsion at contact and weakens the attractive well at intermediate interparticle distances, completely eliminating the latter at high polydispersity index. The reduction in contact repulsion was attributed to polydispersity relieving monomer crowding near the particle surface, especially at high grafting densities. The elimination of the midrange attractive well could be attributed to the longer grafts in the polydisperse graft length distribution that in turn introduced longer range steric repulsion and altered the wetting of the grafted layer by matrix chains. That work demonstrated that at high grafting densities, polydispersity in graft length can be used to stabilize dispersions of grafted nanoparticles in a polymer matrix at conditions where monodisperse brushes would cause aggregation.

Polymer Entanglements
One of the fundamental concepts in the molecular description of structure-property relations of polymer melts is chain entanglement. As the molecular weight of the molecules in a polymer melt is increased, the spatial domain spanned by any given chain increasingly overlaps with those occupied by its neighbors. When macromolecules interpentrate, the term entanglements intends to describe the topological interactions resulting from the uncrossability of chains. The fact that two polymer chains cannot go though each other in the course of their motion changes their dynamical behavior dramatically, without altering their equilibrium properties. Entanglements play a key role in the viscoelastic properties of polymers, as evidenced, for example, by the emergence of a plateau region in measurements of the storage modulus as a function of frequency.
Molecular simulations have confirmed that the overall motion of the chains in a polymer melt is restricted to diffusion along their "primitive paths", which represent the diffusive paths that linear chain molecules follow between their two ends as a result of topological constraints [149]. The advent of computational algorithms enabled direct observation of entanglements that arise in polymeric melts [276][277][278]. Anogiannakis et al. [38] have examined microscopically at what level topological constraints can be described as a collective entanglement effect, as in tube model theories, or as certain pairwise uncrossability interactions, as in slip-link models. They employed a novel methodology, which analyzes entanglement constraints into a complete set of pairwise interactions (links), characterized by a spectrum of confinement strengths. As a measure of the entanglement strength, these authors used the fraction of time for which the links are active. The confinement was found to be mainly imposed by the strongest links. The weak, trapped, uncrossability interactions cannot contribute to the low frequency modulus of an elastomer, or the plateau modulus of a melt.

Insight Obtained from Simulations
Riggleman et al. [279] have carried out a detailed examination of entanglements in a nanocomposite glass. They have conducted Molecular Dynamics simulations of an ideal bead-spring polymer [149] nanocomposite model in which the nanoparticles were dispersed throughout the polymeric matrix. After equilibration in the melt state, all configurations were cooled below their glass transition temperature, where they were subsequently aged using MD for a short period. Finally, a simulation of the creep response of each sample was performed, where tensile and compressive stresses were applied to the glassy specimens. In order to reduce the chains to their primitive paths, these authors employed the CReTA algorithm of Tzoumanekas and Theodorou [278]. During the reduction process, the diameter of the particles is reduced to facilitate slippage of entangled chains past each other, up to the point that further decrease in the diameter of the particles no longer has an appreciable effect on primitive path statistics. The nanoparticles are necessarily frozen in space as the algorithm proceeds.
By examining all particles, one can calculate the distribution of the number of primitive path contacts per particle in the system, shown in Fig. 29(a). The majority of the particles trap at least one primitive path. The entanglements due to polymer chains crossing each other are expected to exhibit little (if any) change during deformation. The only mechanism for the entanglements to disappear is through chain ends slipping past an entanglement junction; such effects are anticipated to be minimal in the glassy state. However, Fig. 29(a) reveals appreciable changes in the distribution of the number of primitive path contacts per particle; the number of contacts per particle increases significantly upon deformation. Figure 29(b) shows how the average number of contacts per particle increases with time for both tensile and compressive deformations. An intriguing, overall physical picture emerges from the heterogenous nonaffine displacements and the particle-induced nucleation of entanglements ( Fig. 29(b)). Figure 29(c) provides the number of entanglements for three particles. The particle that exhibits the largest nonaffine displacements begins with two primitive path contacts, and as the deformation proceeds loses its primitive path contacts. Since that particle was not hindered by any primitive paths, it was able to move throughout the system more easily. In contrast, the particle with the smallest nonaffine displacements (plotted using up triangles) experienced two or more primitive path contacts during the entire deformation. Those entanglements served to trap the particle and forced it to move along them, in an affine manner. Nanoparticles were found to serve as entanglement attractors, particularly at large deformations, altering the topological constraint network that arises in the composite material.
Hoy and Grest [280] performed primitive path analysis [276] of polymer brushes embedded in long-chain melts. All simulations were for a coarse-grained model [149] in which monomers were represented by beads (of size ) connected by springs. The systems studied consisted of long grafted chains of length N = 501 beads, whereas the entanglement length in a melt is approximately N e = 70 [281]. The polymeric matrix studied consisted of melt chains of length P = 1000 beads. As expected, the brush-brush entanglement density, bb e (z), increases rapidly with the grafting density for overlapping brushes. The brush-melt entanglement density, bm e (z), increases also with the grafting density, but even at low grafting densities there is considerable brush-melt entanglement. Moreover, there is clear crossover from dominance of brush-melt entanglements to brush-brush entanglements as coverage increases. Figure 30 depicts brush-brush, brushmelt and melt-melt entanglement densities for three different grafted densities 0.008, 0.03 and 0.07 (in units of −2 ). The peak of bm e (z) is always at z ≃ 15 , but the width of the first peak increases dramatically with increasing grafting density. At low z, the crossover between a preponderance of brush-melt entanglements and preponderance of brush-brush entanglements clearly occurs at 0.03 −2 grafting density. At this coverage, the peaks of the brush-brush and brushmelt entanglement densities are of nearly equal height. For higher coverages, the peak of the brush-brush entanglement density is higher, the reverse of the situation for lower coverages. Summarizing, when surrounded by melt, the brushes entangle predominantly with the melt at low coverage and with themselves at high coverage. The peak of the brush-melt entanglement density is highest at an intermediate coverage, but the integrated areal brush-melt entanglement density continues to increase with coverage for the studied systems.

Experimental Findings
Nanoparticles have been shown to influence mechanical properties, as well as transport properties, such as viscosity. The total number of primitive path contacts per nanoparticle as the system deforms in tension (solid line) or compression (dashed line). c Number of primitive path contacts plotted against the instantaneous strain for three chosen particles as the nanocomposite system deforms. The particle that exhibited the largest nonaffine displacements is represented by left triangles while the one with the smallest nonaffine displacements is plotted using up triangles. (Reprinted from [279], with the permission of AIP Publishing.) Until recently, the commonly held opinion was that particle addition to liquids, including polymeric liquids, produces an increase in viscosity, as predicted by Einstein a century ago [282,283]. However, it was recently found by Mackay and coworkers [208,284,285] that the viscosity of polystyrene melts blended with crosslinked polystyrene particles (and later also with fullerenes and other particles) decreases and scales with the change in free volume (due to introduction of athermal excluded volume regions in the melt) and not with the decrease in entanglement. Later, [285] fullerenes and magnetite particles were found to produce the same non-Einstein viscosity decrease effect.
Micron-sized spherical fillers increase the viscosity of a pure polymer melt from p to a value of predicted by the Einstein-Batchelor law: where is the particle volume fraction [201,219,286]. However, for nanosized fillers, can be reduced or increased relative to the pure polymer [208,284,[287][288][289][290][291][292][293]. While there have been extensive simulations on nanocomposites, a few of them have focused on the importance of nanoparticle addition on flow behavior [290,294].

Insight Obtained from Simulations
Kalathi et al. [196] employed nonequilibrium Molecular Dynamics simulations in order to find out whether the shear viscosity of a polymer melt can be significantly reduced when filled with small energetically neutral nanoparticles. That proved to be the case, apparently independently of the polymer's chain length. Analogous to solvent molecules, small nanoparticles seem to act as plasticizers and reduce the viscosity of a polymer melt. Their simulations allowed them to organize the viscosity data of filled polymer melts as a function of the dimensions of the matrix chains and the particles. Figure 31 Square symbols correspond to ∕ p < 1, diamonds to ∕ p > 1, circles to ∕ p ≃ 1 at low nanoparticle loading, and triangles to the case where an initial increase of viscosity with nanoparticle loading is followed by a decrease. a Experimental data for athermal systems are from [208,284]. Systems above the solid orange line should be miscible. The black "viscosity" line is extrapolated from the simulation findings. b Corresponding plot for dissimilar mixtures. Only the viscosity line is shown, Data are from [287][288][289][290][291][292][293] (Reprinted figure with permission from [196]. Copyright 2001 by the American Physical Society) melts have the same chemical structure [208,284], while chemically dissimilar mixtures are considered in Fig. 31(b) [287][288][289][290][291][292][293]. In Fig. 31(a), Kalathi et al. have also included the miscibility line from ref. [150] suggesting that the experiments correspond to miscible nanoparticle-polymer mixtures. The "viscosity" line drawn from the simulations separates regions where the nanocomposite's viscosity is smaller from those where viscosity is larger than that of the pure melt. For short chains, the viscosity crossover occurs when the nanoparticle size is comparable to R g . In contrast, the limited data for large R g suggest that the line is nearly vertical.
Stephanou et al. [295] introduced a continuum model for polymer melts filled with nanoparticles capable of describing in a unified way their microstructure, phase behavior, and rheology in both the linear and nonlinear regimes. That model was based on the Hamiltonian formulation of transport phenomena for fluids with a complex microstructure with the final dynamical equations derived by means of a generalized (Poisson plus dissipative) bracket. The model describes the polymer nanocomposite melt at a mesoscopic level by using three fields (state variables): a vectorial (the momentum density) and two tensorial ones (the conformation tensor for polymer chains and the orientation tensor for nanoparticles). A key ingredient of the model is the expression for the Helmholtz energy, A, of the polymer nanocomposite. Beyond equilibrium, A contains additional terms that account for the coupling between microstructure and flow. In the absence of chain elasticity, the proposed evolution equations capture known models for the hydrodynamics of a Newtonian suspension of particles. Figure 32 presents the relative viscosity predicted by the model of Stephanou et al. for an unentangled PEO melt with molecular weight M = 1000 g/mol filled with silica nanoparticles of diameter D = 43 nm [296]. Due to the large nanoparticle volume fractions covered in the measurements (up to 50%), neither the Einstein equation nor the Einstein-Berthelot-Green one are applicable [296]. A better choice is the empirical equation proposed by Krieger and Dougherty [297] for dense Newtonian suspensions. It can also be observed that for ≥ 0.27 the data in Fig. 32 exhibit a plateau in the limit of infinitely high shear rates. At those shear rates, flow is so fast that thermal motion cannot destroy the imposed structure (fully aligned molecules) on polymer chains.

Moduli of PNCs
Three different models [298] have been proposed by Einstein [282,283], Eilers [299], and Guth [300] for estimating the enhancement of the shear modulus of composites incorporating spherical particles: with being the volume fraction of particles dispersed and G and G p the shear moduli of the pure polymer and the composite, respectively. Einstein derived his model for small volume fractions of particles, where the enhancement in the shear modulus (or viscosity increase) can be estimated by a linear superposition of the shear distortions arising from individual particles; though this relationship was originally derived for shear viscosity of particle suspensions, it is also applicable to a host of other properties, including the shear modulus of composites. Later, Guth extended this model to higher by accounting for additional shear distortion arising from the interactions between the distortions arising from neighboring particles. Eilers made empirical corrections to Einstein model to account for the dramatic rise in the viscosity of suspensions observed when the volume fraction approaches the closepacking sphere density limit.
Surve et al. [312] employed a combination of polymer mean field theory and Monte Carlo simulations to study the polymer-bridged gelation, clustering behavior, and elastic moduli of polymer-nanoparticle mixtures. Polymer selfconsistent field theory was first numerically implemented in order to quantify both the polymer induced interaction potentials and the conformational statistics of polymer chains between two spherical particles. Subsequently, the formation and structure of polymer-bridged nanoparticle gels were examined using Monte Carlo simulations. These authors used the number distribution of bridges, obtained from their simulations, to quantify the elastic properties of the polymer nanocomposites in the postgel regime. Similar to classical network theories, they assumed that the only contribution to the elastic response of the system comes from the backbone of the percolated network and that the "sol" fraction and the dangling ends of the network do not impact elasticity to the system. They defined the backbone of the percolated network as the percolated cluster, excluding the dangling tails and dangling loops, that can be identified as the largest biconnected component of a percolated cluster. Since, for the case of bridging induced percolation, the interparticle bridges served as the stress bearing bonds between the particles, the enhancement in the elastic modulus was assumed to be proportional to the number of such bridges, at a given volume fraction of particles. Figure 33(a) displays the elastic moduli scaled by a constant factor as a function of particle volume fraction, expressed as − c , with and c being the volume fraction and the percolation volume fraction of the particles, respectively. As observed from the figure, the elastic moduli follow a universal power law scaling, G ∝ − c , with ≃ 1.79. If energetic contributions to elasticity are taken into account, higher elastic exponents appear from ∼2.1 to 3.8 depending on the relative influence of stretching entropy and bending energy [313,314] (Fig. 33(b)).
McEwan et al. [315] predicted the storage modulus by employing a Zwanzig-Mountain relation and Monte Carlo simulations. In parallel, these authors measured the modulus from rheology experiments on samples well characterized with ultra-small angle X-ray scattering. These authors connected particle microstructure to the storage modulus at infinite frequency, G � ∞ , through the Zwanzig-Mountain equation for isotropic molecular fluids [316]: where the first term within the parentheses represents an ideal contribution to the modulus due to the presence of the particles and the second term accounts for the contribution (70) . 33 a Master curve for the scaled elastic modulus obtained from simulations as a function of the particle volume fraction − c . b Scaled elastic modulus for experimental polymer-particle sys-tems. The sources of experimental data are listed in refs [301][302][303][304][305][306][307][308][309][310][311].
(Reprinted from [312], with the permission of AIP Publishing.) due to their interaction and microstructure. The radial distribution of nanoparticles, g(r), was obtained from Monte Carlo simulations of the particles obeying a Mewis-Russel [317] potential for polymer-grafted spheres. Their results are presented in Fig. 34. Silica particles of radii R ≃ 100 nm and 600 nm were synthesized, grafted with hydroxylterminated polydimethylosiloxane (PDMS) chains and finally dispersed in PDMS matrices at volume fractions, c , ranging from 0.02 to 0.65. The experimental measurements are presented alongside the theoretical results in Fig. 34. It can be clearly seen that the storage modulus increases upon the addition of particles, following an almost universal scaling with the volume fraction of the particles. In all cases, the Zwanzing-Mountain predictions are very close to the experiments and to the predictions of the Hall equation of state for solids [318] (at high volume fractions where the materials behave in a solid-like manner).
To provide insights into how polymer-grafted nanoparticles (NPs) enhance the viscoelastic properties of polymers, Hattemer and Arya [319] computed the frequencydependent storage and loss moduli of coarse-grained models of polymer nanocomposites by employing Molecular Dynamics simulations. Figures 35(a) and (b) present the computed G � ∕G � 0 and G �� ∕G �� 0 ratios plotted against for the six polymer nanocomposites containing bare and grafted nanoparticles at low-and high-frequency regimes along with the the moduli predicted by using the theoretical models of Einstein, Eilers and Guth. It is evident that both G � ∕G � 0 and G �� ∕G �� 0 increase with increasing volume fraction, consistent with the trend obtained from the strain distortion models, though the different models differ somewhat from each other. Moreover, the moduli ratio computed from simulations at high frequencies are of comparable magnitude to those predicted by the models, especially the model proposed by Guth, whereas the ratios computed at low frequencies tend to exceed all model predictions. At low frequencies, the computed G � ∕G � 0 ratios are more strongly affected compared to G �� ∕G �� 0 , which is consistent with the expectation that G ′ is more strongly affected by changes in the relaxation times of the polymer chains (quadratic dependence with the Rouse time) as compared to G ′′ (linear dependence).

Fig. 35
Ratio of storage (a) and loss modulus (b) of the three bareand the three grafted-nanoparticle systems to that of pure polymer plotted as a function of the effective nanoparticle volume fraction. The moduli ratios obtained from simulations at low and high frequencies are shown as blue circles and red squares, respectively. Open symbols represent bare nanoparticles, and filled symbols represent grafted nanoparticles. The black solid, dotted, and dashed lines represent predictions from the models of Einstein [282,283], Eilers [299], and Guth [300], respectively. (Color figure online) (Reprinted with permission from Ref. [319]. Copyright (2015) American Chemical Society.) Fig. 36 a Distribution of local C 44 with respect to the distance r from the center of a nanoparticle for three different types of interaction considered at temperature T = 0.5. b Distribution of the local C 44 with respect to the distance r from the center of the nanoparticle for the attractive particle in the melt and glass regime. (Reprinted figure with permission from [320]. Copyright 2005 by the American Physical Society)

Local Moduli
It is now generally accepted that a nanoparticle will perturb the conformation of the polymer around it. However, the question still remains open whether such conformational changes are directly responsible for the mechanical behavior of the polymer, i.e. whether there are any reinforcement or weakening effects of a polymer by a nanometer-sized particles, whether such effects are localized, and if so, what is the extent and the magnitude of that localization. Papakonstantopoulos et al. [320,321] have developed a formalism and applied it to calculate the local mechanical properties of a nanocomposite system in detail. Their coarse-grained, bead-spring Monte Carlo simulations revealed that a glassy layer is formed in the vicinity of the attractive filler, contributing to the increased stiffness of the composite material. Following Yoshimoto et al. [322], the local mechanical properties of the system were determined by discretizing the simulation box into small cubic elements and measuring the internal stress fluctuations within each cubic subdomain [323]. Figure 36(a) shows the local shear modulus as a function of the distance from the surface of the filler. An increase of the local C 44 is observed for the attractive systems in the vicinity of the particle. This pronounced increase may be indicative of the existence of a glassy layer around the particles, even at temperatures above the glass transition temperature (T∕T g = 1.16), which was hypothesized by Berriot et al. [324]. The results of the neutral and repulsive system are more intriguing. The nanoparticle is surrounded by a region of negative modulus which is followed by a second region of higher than the bulk modulus. Figure 36(b) shows the local shear modulus as a function of temperature for the attractive particle. It can be seen that, as the temperature decreases, the shear modulus of the solid-like layer around the particle increases. The thickness of that glassy layer, which is comparable to the radius of gyration of the polymer, also increases. In all cases, far from the particle, the shear modulus decays to the value corresponding to the pure polymer at the given temperature, as expected. Summarizing, it seems that, even above the glass transition temperature, nanoparticles induce the formation of a solid-like layer, whose existence has been invoked to explain experimentally observed increases of the storage modulus in nanocomposites.

Deformation Simulations
Riggleman et al. [325] have examined the response of a polymer and a polymer nanocomposite glass to creep and constant strain rate deformations using Monte Carlo and Molecular Dynamics simulations. These authors found that nanoparticles stiffened the polymer glass, as evidenced by an increase in the initial slope of the stress-strain curve and a suppression of the creep response. Figure 37 shows the stress-strain curves obtained by Riggleman et al. [325] for both the neat and the nanocomposite polymer for both tension and compression at two different strain rates. All curves exhibit similar features: it can be discerned an initial elastic response followed by yield and strain softening when the strain ≃ 0.05. For strains beyond ≃ 0.10 the stress rises again as strain hardening begins. The Young's modulus, E, was obtained by fitting the linear part of the elastic response ( ≤ 0.02) = E . Both under tension and compression, the nanocomposite system was found to be stiffer. Moreover, these authors reported that constant strain rate and constant stress deformations had different effect on the material's position on its energy landscape, in a way that neither the stress nor the strain rate were uniquely indicative of the relaxation times in the material [326].
Chao and Riggleman [327] studied the effect of nanoparticle curvature and grafting density on the mechanical properties of polymer nanocomposites. In their study, they developed a coarse-grained model of a polymer glass containing grafted nanoparticles and examined the resulting effects on the elastic constants, strain hardening modulus, as well as the mobility of the polymer segments during deformation. They found that the elastic constants and yield properties were enhanced nearly uniformly for all nanocomposite systems studied, while the strain hardening modulus depended weakly on the grafted density and the nanoparticle size. Figure 38 shows the mechanical response of the systems studied under compressive deformation at a constant rate, where the measured stress is plotted against the ideal rubber elasticity factor, g( ) = 1∕ − 2 , with being the macroscopic stretch imposed on the specimens. Early in the deformation, the polymer nanocomposites exhibit an elastic response (g( ) < 0.15) followed by yielding and strain softening. Finally, at larger stretches (g( ) > 0.3), the polymer glasses enter the strain hardening regime, and the stress resumes an increasing trend as strain continues to grow. These authors decomposed the stress calculated in their simulations into its components. The normalized contribution from the non-bonded interactions between the nanoparticles and the grafted chains is presented in Fig. 38(a). It can be seen that the interaction between the nanoparticle and the grafted chains increases with the particle size, and its effect is only observed in the strain hardening region. This finding is coherent with the expected depletion of the matrix from particle surfaces with increasing particle size. Similarly, the stress contributions from the non-bonded interactions among the grafted chains is presented in Fig. 38(b). For low grafting densities (0.05 and 0.1) the non-bonded interaction between beads belonging to grafted chains does not contribute significantly to the stress increase during strain hardening. However, for higher grafting densities (0.2 and 0.4), the nonbonded interactions between the grafted chains contribute significantly to strain hardening. In contrast to the obvious dependence of the non-bonded component of the stress tensor to the grafting density and particle radius, the bonded component of the stress tensor did not exhibit significant changes in behavior across the various nanocomposite systems investigated by Chao and Riggleman [327].
Hagita et al. [328] performed coarse-grained Molecular Dynamics simulations of nanocomposite rubbers with spherical nanoparticles on the basis of the Kremer-Grest [149] model. Figure 39 shows the stress-strain relations of a small mesh cross-linked polymer network for three nanoparticle-polymer interactions (repulsive, slightly attractive and attractive). There are clear differences between the repulsive and the attractive cases. However, both cases with attractive nanoparticle-polymer interactions seem to behave similarly, probably due to the few contacts existing between the nanoparticles. Thus, the effect of the exact interaction strength on the stress-strain relations is minor. When the nanoparticles are trapped and fixed in a cross-linked polymer network, the number of contacts between them is expected to increase for a larger elongation ratio due to the compression in the directions perpendicular to the elongation axis. These authors [328] have also calculated the two-dimensional scattering patterns of nanoparticles during the elongation of the network. For strain levels >50% they observed a spot pattern in the structure factor and a two-point bar pattern in the scattering intensity.

Concluding Remarks
We have presented a detailed, result-driven, review of research to address the fundamental problem of PNCs by implementing computer simulations at different levels of description. It is unfeasible through the use of a single simulation technique to capture all the relevant physics of the problem. On the one hand, fully atomistic Molecular Dynamics (MD) can account for the chemical interaction between nanoparticles and the polymer matrix. However, due to the computational demands of atomistic MD, the polymer matrix has to be comprised of oligomers rather than entangled polymers. On the other hand, coarsegrained methods can provide us with an understanding of the underlying phenomena. However, due to the lack of explicit chemical information, coarse-grained methods should be carefully parameterized based on findings of more detailed simulations. In any case, even the wide spectrum of molecular simulation methods developed till now cannot fully capture the macroscopic behavior of PNCs. Coupling molecular simulations to continuum calculations is the way to achieve this [329][330][331].