Elasticity of polymeric nanocolloidal particles

Softness is an essential mechanical feature of macromolecular particles such as polymer-grafted nanocolloids, polyelectrolyte networks, cross-linked microgels as well as block copolymer and dendrimer micelles. Elasticity of individual particles directly controls their swelling, wetting, and adsorption behaviour, their aggregation and self-assembly as well as structural and rheological properties of suspensions. Here we use numerical simulations and self-consistent field theory to study the deformation behaviour of a single spherical polymer brush upon diametral compression. We observe a universal response, which is rationalised using scaling arguments and interpreted in terms of two coarse-grained models. At small and intermediate compressions the deformation can be accurately reproduced by modelling the brush as a liquid drop, whereas at large compressions the brush behaves as a soft ball. Applicable far beyond the pairwise-additive small-strain regime, the models may be used to describe microelasticity of nanocolloids in severe confinement including dense disordered and crystalline phases.

deformation response virtually independent on the number and length of chains in the SPB. Our predictions can be explored experimentally either directly in force-deformation measurements of SPBs 19 or indirectly by analysing the structure of their condensed phases and aggregates 20 .

Results
We use molecular dynamics (MD) simulations to study the deformation of an SPB consisting of a small hard colloidal particle grafted with a polymer brush of f linear-chain arms each consisting of N c monomers; temperature was fixed at  = / T k B where  is the energy scale of the monomer-monomer repulsion (see Methods). The SPB is confined to a slit formed by parallel immobile walls (Fig. 1a) and its deformation is quantified by the ratio of central lateral extension and indentation denoted by ζ. To evaluate it, we computed the dimensions of the SPB based on monomer densities projected on the axes of the coordinate system centered at the SPB and oriented such that the z axis is perpendicular to the walls; note that the projected densities , , c c x y and c z are measured in units of 1/length rather than in units of 1/volume as the usual monomer density c. SPB half-thickness is given by the transverse semiaxis evaluated in absence of walls. We find that our ζ is more meaningful than its analogue based on the eigenvalues of the radius of gyration tensor, which carries a larger numerical error at small compressions. Figure 1b shows the reduced central lateral extension ζ for thin (N c = 30) and thick (N c = 50) SPBs with functionalities f = 30, 40, 50, and 60; the screening length κ −1 of the Yukawa wall potential is a small fraction of the radius of gyration of an isolated SPB R g 0 ( . R 0 05 g 0 in the thin and . R 0 04 g 0 in the thick SPB). The eight datasets plotted against reduced slit width / L R 2 g 0 collapse surprisingly well despite considerable variation of chain length N c and functionality f, and they reveal three distinct deformation regimes.
At small compressions the MD results are rather scattered due to smallness of both numerator and denominator in equation (1) but they still reveal a knee-like increase of ζ clearly visible in the two f = 60 datasets replotted in Fig. 2a. This behaviour is additionally confirmed by the self-consistent field theory (SCF; see Methods). SCF results shown in the inset to Fig. 1b agree well with the MD data but are considerably less noisy, emphasising the knee-like onset of deformation.
The narrow small-compression regime is followed by a broad, virtually linear variation of ζ extending from reduced slit width / L R 2 g 0 of about 1.6 down to about 0.7. At / ≅ . L R 2 05 g 0 , the SPBs undergo a transition to the large-compression regime characterised by a steep increase of ζ upon compression.  (2)] is plotted with the dashed line faded at / L R 2 g 0  1.5 to emphasize that it is valid only in narrow slits. Inset: ζ obtained from SCF theory for the same f and N c plotted using the same color code; also included is equation (2)  The effective SPB diameter can be defined by the onset of deformation where the slit serves as a vernier caliper. The f = 60 MD data in Fig. 2a  ; this is the effective diameter. Since a linear polymer is more anisometric than an SPB, it must have a larger ratio / ⁎ D R g 0 so that our estimates are reasonable. Finally, the auxiliary Yukawa potential used in simulations introduces an effective slit width smaller than the nominal L by about κ / 2 . Thus the true diameters of the (N c = 30) and the (N c = 50) SPB are smaller than the above D * by about . R 0 10 g 0 and . R 0 08 g 0 , respectively. The remarkable universality of the SPB deformation must stem from a basic feature of polymers, and it is instructive to begin understanding it by a scaling-theory estimate of ζ for a single linear chain. In severe confinement 22 , the in-plane diameter of the chain is ∼ − / / D L N 1 4 3 4 whereas its transverse size is ∼ ⊥ D L, which gives irrespective of chain length. The dashed line in Fig. 1b shows that this result is in good semiquantitative agreement with the MD data although it relies on the blob picture of a linear chain with excluded-volume interactions rather than on the geometrically more involved blob analysis of the deformation of a multiarm brush 23 . The agreement suggests that it is worthwhile seeking a coarse-grained interpretation that does not depend on the details of the macromolecular architecture, and here we propose two complementary continuum theories.
Here γ is the surface tension and V 0 is the reference volume where the pressure within the drop given by the Murnaghan equation of state 24 vanishes in absence of surface tension so that χ T is the isothermal compressibility of the drop at p = 0 and γ = 0. First proposed for hydrostatic compression of elemental substances 24 and previously considered in models of compressible rubber 25 , this equation captures the essential physics of simple fluids. The ideal-gas-like term ensures that the pressure diverges at small volumes, thereby phenomenologically accounting for the repulsive interparticle forces, whereas the negative, volume-independent term represents the effect of the cohesive interactions.
Equation (3) contains a single surface term although it could be split into two contributions corresponding to the free surface of the brush and to the brush-wall contact zone. Such a generalization is justified even in inert walls studied here, yet we stick to the more transparent single-surface-tension variant of the model because it already offers a very accurate interpretation of the shape of the confined SPB as shown below. Within this model, the deformation of the drop is controlled by a single dimensionless parameter T , the term in the parenthesis representing the Laplace pressure. Unlike in ordinary liquids, in which the drop size R * can be arbitrarily large, for spherical polymer brushes R * is a length-scale determined by the brush architecture, f and N c . Accordingly, it is physically meaningful to treat Ψ rather than  EW as an intrinsic material parameter, and we refer to it as the reduced Egelstaff-Widom length. It is also convenient to define the deformation free energy, F LD def , as the difference between the liquid-drop free energy of a confined drop and that of a free drop of radius R * . From equation (3) we obtain where V * and A * are the volume and surface area of the drop of radius R * , respectively, and the energy scale U is given by Note that in this way, the reference lengthscale R 0 has completely dropped out as it should. The model is fully described by two parameters, the reduced Egelstaff-Widom length Ψ which sets the shape of the deformation and U which sets the overall scale of the energy penalty to compress the drop.
In Fig. 2a we compare the f = 60 MD data to the Ψ = 0.6 liquid-drop ζ calculated from the semiaxes of the drop which were obtained numerically using Surface Evolver 27 (see Methods). The model reproduces very well the small-and the intermediate-compression regimes at slit widths ranging from 100% down to about 30% of the effective diameter ⁎ D . Note that Ψ is essentially the sole fitting parameter as the MD data are consistent only with very limited variations of ⁎ D ; here we used = .
In the large-compression regime the model is no longer suitable. This is shown by the continuation of the best-fit Ψ = .
0 6 liquid-drop ζ at small slit widths / < . L R 2 05 g 0 plotted in Fig. 2b, which underestimates waist extension increasingly more as L is decreased and leads us to the condition that in this regime the behaviour of the SPB is qualitatively different. Soft-ball model. To reproduce the deformation at large compressions, we turn to an alternative model where the SPB is represented by an elastic solid sphere. The underlying rationale is that in the large-compression regime, chain fluctuations are reduced considerably by a combination of topological restrictions due to grafting and confinement, implying that the brush may well behave effectively as a solid rather than as a liquid.
Among the several types of elasticity of isotropic media, we choose the modified neo-Hookean theory originating in the statistical thermodynamics of a three-dimensional polymer network 28 . The corresponding free energy density reads 29 where Y is the Young modulus and ν is the Poisson ratio whereas I 1 is the first invariant of the isochoric part of the Green deformation tensor F F T , F being the deformation gradient, and = J F det . Like in the liquid-drop model, it would be reasonable to complement the soft-ball elastic energy by the surface energy. Yet within the minimal framework this is not needed because shear elasticity alone resists shape deformation. To keep the model as simple as possible, we opt to not include the surface tension so that the ball shape is controlled solely by the Poisson ratio ν, whereas the deformation energy scales as ⁎ YD 3 . The equilibrium shape of the ball at a given L was found by minimizing the elastic energy using the FreeFEM+ + package 30 (see Methods). We varied the Poisson ratio so as to obtain an optimal agreement at large compressions, finding that the best-fit value of ν is 0.3 ( , can be understood by visualising a slightly compressed sponge ball where stresses due to compression are localized right at the two walls and do not reach into the bulk, and the waist diameter is thus barely increased. In a liquid drop, on the other hand, any increase of pressure upon compression is communicated across all of the volume and thus the waist extension relative to compression is considerable even at small compressions.
Density profiles. Our theoretical description of the MD shape deformation data combines the predictions of the two models, the liquid-drop/soft-ball transition being at / ≅ . L R 2 05 g 0 . A more detailed insight into the transition as well as the SPB structure itself is provided by the monomer density profiles in Fig. 3  g 0 and 1.9 for the MD and the SCF results, respectively. At this slit width, the agreement of the MD and the SCF density profiles is rather good as expected 31 , the only difference being the slightly larger distance between the 1% SCF density isoline and the theoretical SPB surface which can be attributed to the fact that the effective diameter predicted by SCF is larger than that obtained by MD. The location of the transition from the inner, dense region of the SPB and the outer dilute shell with an approximately exponential density profile coincides with the theoretical surface; on the linear scale of Fig. 3, the small-magnitude exponential tail is seen as an almost straight vertical segment of the red curves. Also visible is the effect of the auxiliary Yukawa potential used in MD simulations which produces a depleted subsurface layer close to either wall, and a few artifacts of the cubic lattice used in SCF theory (say in the kidney-shaped 50% isoline).
As the SPB is compressed, the agreement of the MD and the SCF results is gradually poorer. In part, this is due to the approximate treatment of the excluded-volume interactions in the SCF theory 31 and in part it can be related to the slightly different effective SPB diameters and to the different wall types used. This is best seen in the 20% slit-width case where the SCF midplane radial density profile is representative of the density in the whole SPB whereas the MD profile varies considerably with the distance from the midplane. Still the MD midplane profiles show that the boundary of the inner dense region of the SPB correlates with the theoretical soft-ball/liquid-drop surface. Also notable is the development of the shoulder-like density profile in very compressed SPBs signaling a qualitatively different behaviour compared to small and intermediate compressions, which is consistent with the soft-ball/liquid-drop transition.
Deformation energy. The coarse-grained picture is completed by comparing the MD deformation energy to those of the two models. To this end, the effect of the Yukawa wall potential is taken into account by recognizing that at small slit widths  L R 2 g 0 the Yukawa potential penetrates across the whole slit. The corresponding energy increase can be estimated by the magnitude of the potential in the center which is proportional to Figure 4 shows the MD energy for the N c = 30 and the N c = 50 f = 60 SPBs as well as the fits obtained by combining the liquid-drop and the soft-ball energies with the above Yukawa term. The agreement is very good over almost three orders of magnitude, the only systematic but insignificant deviation being the behaviour at > ≅ .
⁎ L D R 3 42 g 0 where the liquid-drop and the soft-ball energies vanish whereas the MD energy remains finite. A close inspection reveals a minute discontinuity in the deformation energy at the liquid-drop/soft-ball transition, which may be due to the approximate treatment of the Yukawa potential.
The best fits of deformation free energy in Fig. 4 fix the characteristic energy scales of the two models, U [equation (7)  Scaling theory. Our numerical analysis revealed a remarkable collapse of SPB deformation data. The eight datasets shown in Fig. 1b cover a broad range of experimentally relevant SPB functionalities and two chain lengths large enough so as to ensure that the reported behaviour is not affected by the monomer size. Here we provide a scaling-theory interpretation of the data collapse, which suggests that our observations are universal.
Our starting point is the des Cloizeaux formula for the osmotic pressure Π of semidilute polymer Ta   B 15 4 9 4 , where φ is the monomer concentration and a is the monomer size 32 . This result can be used to calculate the reduced Egelstaff-Widom length [equation (4)] recast as χ χ Ψ = λ Π ≅ Π Ψ T T so as to emphasise that the SPB surface tension is related to the osmotic pressure by the Young-Laplace equation; here we note that for Ψ = . and This implies that χ ∼ /Π 1 T and thus the reduced Egelstaff-Widom length is a quantity of order unity that does not depend on functionality nor on chain length, This finding is in very good agreement with the data in Fig. 1b and it   chain length. Note that this prediction is valid for any power-law equation of state rather than just for the des Cloizeaux formula, and it can be attributed to the self-similar chain structure. The small deviations from the perfect data collapse may be due to the presence of the colloid because its size is the same in all SPBs studied here whereas the radius of SPB diameter varies with N c and f. Within the same theory, we can further extract the dependence of the energy scale U [equation (7)] on the brush parameters. Since This is a very rewarding result, as it coincides with the scaling of both the pairwise effective interaction potential between star polymers 33 and between a star and a single wall 34 . Indeed, in the limit of weak compressions, the overall deformation energy penalty paid by the brush should be pairwise additive; accordingly, the elastic theory put forward should reproduce the dependence of the interaction on f, as it does. Moreover, and once more in agreement with the aforementioned effective potential, there is no dependence of the energy scale on N c . This finding is also confirmed in our simulations: as can be seen in Fig. 2b, the deformation energy curves for the two brushes of the same f but different N c run very close to one another for a broad range of deformations in which the liquid drop model is valid. The small deviations seen are of order 10% and can be attributed to the rigid core which makes the thinner brush effectively less compressible.
Finally, we can now provide an independent estimate of the surface tension and compressibility, based on purely theoretical arguments, and compare with the values obtained by the fit. From equations (4,7,11,12), we readily obtain the scaling laws

Discussion
The universal value of the reduced Egelstaff-Widom length Ψ in SPBs is the most important result of the scaling theory obtained by approximating the brush by a star polymer and disregarding the hard colloidal core. Had the core been taken into account, the overall SPB compressibility and thus Ψ should be smaller, leading to a departure from universality. A sizable core is also expected to shift the liquid-drop/soft ball transition to narrower slits by excluding the chains from the centre where a large enough local monomer density could otherwise be reached for geometrical reasons. On the other hand, Ψ can also be controlled by intermonomer interaction and by the solvent, affecting both effective surface tension and compressibility of the SPB. Thus it seems reasonable to theoretically explore the liquid-drop model more closely for a broad range of Ψ covering not only the regime directly applicable to our SPBs but also the incompressible limit at Ψ  1 where the brush is deformed at constant volume and the tension-dominated limit at Ψ  1 where the shape of the SPB is dictated primarily by surface tension. 0 01 drop expands very dramatically, assuming a pronounced pancake-like shape at slit widths below about 40% of D * , whereas the in-plane diameter of the tension-dominated Ψ = 100 drop is only slightly larger than D * . This insight is quantitatively elaborated in Fig. 5b which shows the theoretical reduced central lateral extension ζ of a diametrally compressed liquid drop for Ψ = . , . , , , 0 01 0 1 1 10 and 100. In tension-dominated drops with Ψ  1, ζ is small and barely depends on the ratio of slit width and reference drop diameter / ⁎ L D whereas in incompressible drops of Ψ  1, ζ diverges at vanishing slit widths. In the incompressible limit Ψ = 0, drop shape can be approximated by a jelly-doughnut shape with flat faces pressed against the walls and a rim of semi-circular cross-section, allowing ζ to be calculated analytically. The result (dashed line in Fig. 5b) is characterized by a − / L 1 2 divergence in thin slits as well as a finite ζ at infinitesimally small compressions, which is qualitatively consistent with the numerically obtained knee-like onset of deformation (Fig. 1b). The red shading in Fig. 5b represents the region around Ψ = .
0 6 characteristic of our SPBs, showing that they are halfway between the incompressible and the tension-dominated regime.
The differential response of small-and large-Ψ drops upon compression is even more pronounced in drops confined from all sides rather than just between two walls. This can be illustrated by a drop within a cube-shaped container (Fig. 5c). The nearly incompressible Ψ = .
0 01 drop develops large facets pressed against container walls even at small compressions-as soon as cube edge length L reaches 80% of the reference diameter D * , facets add to 90% of the total drop area. On the other hand, in the tension-dominated Ψ = 100 drop this happens only at / = % ⁎ L D 40 , which corresponds to an eightfold larger density. These differences are important because shape deformation similar to faceting is expected to take place in dense suspensions. Here SPBs are pressed against each other, which too can be described by the liquid-drop model (or, more generally, the combined liquid-drop/soft-ball model) much like repulsion between a star polymer and a hard wall also covers the interaction between two star polymers at small 0 01 and 100 confined to a cubic cavity where the difference between the incompressible and the tension-dominated regimes are readily seen even at small compressions. As the size of the cavity is decreased, the almost incompressible Ψ = .
0 01 drop develops a faceted shape as soon as the cube edge length reaches about 80% of D * whereas the tension-dominated drop maintains a more rounded shape down to cube edge lengths of about 40% of D * . In panels a and c, confining walls are only outlined in the / = % ⁎ L D 100 for clarity.
center-to-center distances 34 . Thus our coarse-grained framework may be used to interpret the structure of crystalline 35 , quasicrystalline 35,36 , and glassy 37,38 nanocolloidal systems as well as certain aspects of their unique rheology 39 , the value of the reduced Egelstaff-Widom length Ψ being the key material parameter.
An interesting direct application of our results addresses the morphologies of aggregated polymer-grafted nanoparticles, which include string-and sheet-like aggregates in addition to the more common spherical bulk-like clusters and dispersed solutions 20 . By using the detailed description of brush deformation in the slit we can readily analyse the lateral repulsive barrier in string-like aggregates, thereby complementing the scaling-theory 40 and patchy-particle 41 interpretations of their stability. Like in the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory, the barrier appears because of competition of the strong but short-ranged van der Waals attraction and a repulsive interaction between them, except that the latter is due to brush deformation rather than electrostatic.
To see how this happens, consider two isolated nearby SPBs (Fig. 6a). Brought into contact by the van der Waals attraction (assumed for simplicity to be dominated by interaction between nanoparticle cores), they will interpenetrate each other as long as the effective elastic modulus of the brush is not too large. Elastic repulsion, which is approximately proportional to a power of indentation with an exponent of 2.1 (inset to Fig. 2b), will merely reduce the magnitude of attraction and the SPBs will eventually form a core-to-core packed dimer with somewhat deformed brushes. In a similar fashion, a third, fourth… SPB may approach the aggregate, and it is intuitively clear that the most favorable point of attachment is at the end of the string where the central SPB is thinnest 40 .
Let us substantiate this expectation by an argument based on the transverse extension of SPBs in the string, which are deformed as if they were confined to a slit. The extension depends on the lengthwise compression dictated by the ratio of core and brush diameters  (Fig. 1b). For another SPB approaching the string from the side, the distance upon contact is thus larger than D * by δ / = % ⁎ D D 2 13 (Fig. 5b). This may not seem much but two important factors come into play that both contribute to the stabilisation of the string. Firstly, the van der Waals interaction is a very steep function of separation and the Hamaker theory predicts that at this larger distance, the core-core attraction between an approaching SPB and a member of the string is only about 46% of the attraction between two SPBs separated by D the confinement of all f chains to a pancake-shaped quasi-two-dimensional region results in an effectively larger functionality, and thus in a decrease of brush compressibility χ ∼ − / f T 3 2 . Provided that χ T is large enough, this may lead to a total potential that is qualitatively different from that in Fig. 5a because the brush-brush repulsion produces a barrier separating the unbound state from the global minimum of core-to-core packed SPBs. If higher than the thermal energy, the barrier prevents the formation of sheetand bulklike aggregates, thereby stabilising strings by a patently many-body mechanism. Alternatively, if the barrier is still not high enough, then aggregation of brushes will also start in a second dimension, leading to the formation of sheets and pushing the SPBs out in the direction perpendicular to the sheets. The two stabilising mechanisms mentioned above will be activated again, potentially preventing the appearance of bulk aggregates. On the other extreme, for  f 1 and very long chains, the chain length itself provides an effective stabilisation against coagulation, and a well-dispersed system results.
In conclusion, the elasticity of complex nanoparticles can be interpreted using coarse-grained notions from continuum mechanics, and the excellent agreement of the simulation results with the liquid-drop model provides a convincing microscopic support for a previously proposed theory of nanocolloidal crystals 42 . Our stripped-down framework can be refined in various ways, say by introducing a differential tension of the free and the contact surface of the brush or by a position-dependent Young modulus which would better correspond to the inhomogeneous monomer density, and it may be extended to non-spherical and polyelectrolyte SPBs 43 which may behave quite differently from our self-avoiding brushes as well as to semiflexible and stiff chains 44 . In addition, the liquid-drop model could be used for the description of viscoelastic effects expected in a dynamic rather than static deformation of the SPB studied here-in this case, one would have to solve the equation of motion for the drop treated as a viscous liquid. Such a generalization may well be relevant for the rheology of suspensions of SPB-like soft colloids, e.g., star polymers 45 .
The most intriguing conceptual conclusion reached is that the deformation behaviour of our SPBs can be encoded solely by the reduced Egelstaff-Widom length Ψ, a dimensionless quantity based on the product of surface tension and compressibility which is known to have a very similar value in many simple liquids 26 . Moreover, the value of Ψ is universal, i.e., independent of the number of chains in the brush and their degree of polymerisation. These findings call for further verification. In turn, it would be interesting to see whether our model also applies to other polymeric nanocolloidal particles such as dendrimer and diblock copolymer micelles and if so, how does the reduced Egelstaff-Widom length depend on the macromolecular architecture.

Molecular dynamics simulations.
Our spherical brush consists of bead-and-spring chains terminally grafted onto a small colloidal particle such that the anchor points are fixed and distributed uniformly across the particle. The main parameters of the SPB are the number of chains f and the number of monomers per chain N c ; colloid diameter is σ = D 8 core LJ where σ LJ is the monomer diameter. Like in a related study 46 , steric monomer-monomer and monomer-colloid interaction is described by the Weeks-Chandler-Andersen (WCA) repulsion 47 and the bonds are modelled by the finite extensible nonlinear elastic (FENE) potential tuned such that the bonds do not cross. Our implicit-solvent scheme and this particular choice of interactions are both consistent with the good solvent conditions studied here. The parameters of the FENE potential used were  σ = / k 30 LJ 2  ( being the strength of the WCA repulsion) and σ = . where m is the monomer mass) and then at any given slit width the SPB was equilibrated for 2 × 10 6 steps followed by 8 × 10 6 measurement steps.
Self-consistent field theory. The self-consistent field (SCF) theory was implemented on a cubic lattice using the Scheutjens-Fleer scheme 46,50 . We chose to represent the walls by a hard rather than Yukawa potential like in MD simulation. This allowed us to analyse the effect of confinement and the SPB geometry as transparently as possible and facilitated comparison with the liquid-drop and soft-ball models.
Continuum mechanics models. In both models, walls were represented by hard constraints so as to study the model in as simple a geometry as possible; also neglected is the SPB core. The equilibrium shape of the SPB represented by the liquid drop was obtained using the Surface Evolver package 27 where the drop is represented by a triangulated surface. The shape was found by relaxation mimicking overdamped motion of the vertices to minimize the combined bulk and surface energy [equation (3)]; the typical number of nodes on the surface was about 3000.
Scientific RepoRts | 5:15854 | DOi: 10.1038/srep15854 The soft-ball model was solved by minimizing equation (8) using the finite-element method implemented within the FreeFEM+ + package 30 . The displacement field was computed on a 3D mesh of about 40000 nodes using a multifrontal LU factorization solver UMFPACK. At any slit width, equilibrium was achieved using Newton-Raphson iteration scheme. To enforce the hard-wall constraints, the magnitude of the parabolic wall potential was gradually increased until the corresponding energy penalty was negligible.