Interaction between Ideal Neutral Nanogels: A Monte Carlo Simulation

Article pubs.acs.org/Macromolecules

Interaction between Ideal Neutral Nanogels: A Monte Carlo Simulation Study Silvia Ahualli,† Alberto Martín-Molina,† José Alberto Maroto-Centeno, and Manuel Quesada-Pérez*,‡ †

Departamento de Física Aplicada, Facultad de Ciencias, Universidad de Granada, 18071 Granada, Spain Departamento de Física, Escuela Politécnica Superior de Linares, Universidad de Jaén, 23700, Linares, Jaén, Spain

‡

S Supporting Information *

ABSTRACT: The interaction between nanogels is a central question in the field of soft matter which has been scarcely studied. In fact, effective potentials for nanogels are less advanced than for other colloidal particles. The soft-sphere potential and the Hertz potential are two theoretical formalisms used until now by some authors to quantify forces between nano- and microgels. Accordingly, in this work we have performed explicit coarsegrained Monte Carlo simulations to find out if these models can capture the interactions between overlapping neutral nanogels. To this end, pairs of nanogels with different number of monomers per chain have been simulated, and the corresponding effective interaction potentials have been calculated as a function of the distance between their respective centers of mass. Then, our simulation results have been used to analyze the functional form of the soft-sphere and Hertz potentials. Both are too simple to describe nanogel interactions when these nanoparticles overlap to a significant extent and can only be employed as a first approach in the regime of small deformations. The simulations performed here also reveal that the interaction can be described within the Hertz model in terms of a bit greater diameter than the geometrical one. In addition, simulations together with theoretical calculations have also performed to deepen into the interaction strength parameter included in the Hertz model. This parameter is a function of the nominal particle diameter as well as Young’s modulus and Poisson’s ratios, and it might exhibit a nonmonotonic behavior that depends on the number of monomers per chain in the nanogels. The results presented here are also analyzed modifying the exponent of the classical Hertz model.

■

INTRODUCTION Nanogels are soft cross-linked polymeric networks with sizes ranging from a few to hundreds of nanometers. In the upper limit, they are also known as microgels. These nanoparticles have been widely studied in the past decades due to their versatility and high sensitivity to external stimuli such as temperature, pH, ionic strength, solvent, external stress, or electric field. These features together with the biocompatibility of some nanogels allow their use in biomedical applications.1 Another interesting property of nanogels is their ability to act as rheological regulators, with potential applications in industry.2 From a theoretical viewpoint, these nanoparticles bridge the gap between classical hard colloids and ultrasoft polymeric colloids, such as star polymers. Nanogels are also employed in experiments addressing condensed matter problems such as structure formation, dynamics, or phase transitions. In any case, our ability to control processes in which nanogels are involved relies on the precise knowledge of the effective interactions between them, which are essential constituents for any attempt to theoretically model and understand these systems. Unfortunately, there are relatively few papers dealing directly with the effective interactions between nanogels.3−6 In © XXXX American Chemical Society

particular, there are few works on the forces between two nanogels when they come into contact, usually referred to as elastic or steric forces. For instance, Richtering and co-workers concluded from different experiments (including rheological studies) that the hard-sphere interaction can be a helpful approximate model for particle volume fractions below 0.3, but the soft-sphere potential is more appropriate for higher concentrations, when interparticle distances are much shorter and microgels can even overlap.7,8 The functional form of the corresponding effective soft-sphere potential is ueff (r ) = εn(σn/r )n

(1)

where r is the separation between the centers of mass of two nanogels, σn plays the role of an effective (or nominal) soft particle diameter, n is the power-law exponent describing the potential stiffness, and εn is the associated interaction strength.3 It should be mentioned that for large n values this potential is almost identical to the hard-sphere potential. Wu et al. Received: October 28, 2016 Revised: February 23, 2017

A

DOI: 10.1021/acs.macromol.6b02333 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 1. (a) Two nonoverlapping nanogels whose centers of mass are at a distance of 26 nm. (b) The same nanogels after approaching (their centers of mass are now at a distance of 20 nm). Green beads symbolize the cross-linkers molecules.

concluded that at relatively low particle volume fractions and below the lower critical solution temperature (LCST) the hardsphere potential can justify the freezing behavior (in agreement with Richtering and co-workers’ conclusions), but an attractive term is required when microgel particles collapse.9 In contrast, recent experiments performed by other authors have shown that the spatial ordering of suspensions of swollen microgels is in very good agreement with a model of particles interacting through the so-called Hertz (or Hertzian) potential, which describes the interactions between elastic spheres.5,6 Its functional form is 2.5 ⎧ ⎪ ε (1 − r / σ ) r < σH H H ueff (r ) = ⎨ ⎪ r ≥ σH ⎩0

However, experiments have shown that this theory is valid for small deformations of a nonadhesive elastic sphere against planes or other particles.11 In particular, it has been successfully applied to soft particles such as vesicles11,12 and microgels.6,13,14 For instance, Schmidt et al. performed force−distance measurements using an atomic force microscope (AFM) to characterize thermoresponsive PNIPAM microgel films.14 Therein, force curves were plotted in log−log scale and fitted with the Hertz model to extract the elastic modulus. However, despite the fact that the Hertz model is acceptable for small deformations, some authors claim that the lack of practicable nonlinear elastic contact models frequently compels the inappropriate use of Hertzian models. More specifically, they would have been improperly used in analyzing phenomena of indentation data of poly(vinyl alcohol) (PVA) hydrogels under major deformations15 or in studying nonlinear rheology at high shear-rate of soft particle glasses.16 In this sense, Burmistrova et al. performed AFM indentation measurements for individual poly(N-isopropylacrylamide-co-acrylic acid) microparticles.17 Data obtained therein were fitted by using a modified form of the Hertz model suggested by Dimitriadis et al.18 to calculate the elastic modulus individually at each indentation depth.

(2)

Here σH is again a nominal particle diameter, and εH is the associated interaction strength. The classical Hertz model has been extensively used to study elastic deformation of spherical solids for more than one century. This model assumes that a homogeneous and isotropic sample undergoing small deformations in which geometric nonlinearities are not considered. These assumptions can be questionable for soft particles.10 B


Article

Macromolecules

(Ntot) and the radius of gyration (Rg), grow with Nm. Such properties are summarized in Table 1.

According to this model, force−indentation curves where there is contact for sure can be fitted with a power law whose exponent is within the interval [1.5, 2.5], being 1.5 the exponent (for force) corresponding to the Hertz model itself.17 In any case, the experiments performed by Richtering, Wu, and Mohanty and their respective co-workers were carried out with PNIPAM-based microgels,6−9,19 but different methods of synthesis were employed. For instance, Senff et al. synthesized their microgels by means of emulsion polymerization in the presence of potassium peroxodisulfate,7 which might lend some ionizable groups and electric charge to the polymer network, whereas Mohanty et al. synthesized their microgels by freeradical precipitation polymerization.6 Consequently, the inner structure and the electric charge of these real polymer networks might not be the same exactly, which complicates a comparative analysis of results. In this sense, computer simulations are extremely helpful in order to research the fundamental properties of ideal and tailormade systems (whose main advantage is the control and the accurate knowledge of their composition and inner structure). In fact, the effective interaction potential has been computationally determined within coarse-grain models for self-avoiding polymer chains,20,21 ring polymers,22−24 polyelectrolyte- and polymer-coated colloids,25,26 or mixtures of self-avoiding star polymers and hard colloids27 (cited just as a few illustrative examples). Accordingly, our main objective is to find out to what extent the soft-sphere potential and the Hertzian potential can properly describe the elastic (steric) interaction between overlapping neutral nanogels. Apart from the choice of the functional form, the effective potential between nanogels might require the identification of an effective diameter (for instance, σH in the Hertz potential) with a measurable quantity of the real system (such as the hydrodynamic diameter), as Heyes et al. have claimed.3 In our case, we might wonder: Can σH be identified with a characteristic length obtained from simulations? Such identification is not a trivial matter for soft particles. Therefore, simulation studies on model systems might be able to help identify the relationship between the effective diameter appearing in analytical expressions and other characteristic lengths describing nanogels.

■

Table 1. Some Properties of the Nanogels Simulated in This Work no. of monomers per chain, Nm

total no. of monomers per nanogel, Ntot

4 8 12 16

466 866 1266 1666

radius of gyration, Rg (nm) 4.62 6.54 8.30 9.67

± ± ± ±

0.10 0.12 0.09 0.06

radius of the nanogel, R (nm) 5.97 8.45 10.71 12.48

± ± ± ±

0.13 0.15 0.12 0.08

Interactions. The short-range repulsion between any pair particles due to excluded volume effects can be modeled by means of purely repulsive Weeks−Chandler−Andersen potential:39−43 ⎧ 6 ⎛ 12 ⎞ ⎪ 4εLJ⎜ d − d + 1 ⎟ r ≤ 12 4⎠ ⎝r r6 uLJ(r ) = ⎨ ⎪ ⎪ ⎩0 r>

6

2d

6

2d

(3)

where r is the center-to-center distance between a given pair of particles, εLJ = 4.11 × 10−21 J, and d is the monomer diameter. The same size was assumed for monomer units and cross-linker molecules (d = 0.7 nm). The interaction connecting monomer units and cross-linkers with their neighbors was modeled by harmonic bonds30,31,44 k bond (r − r0)2 (4) 2 where kbond is the elastic constant (kbond = 0.4 N/m) and r0 is the equilibrium bond length (0.7 nm). Generating the Initial Configuration of Two Interacting Nanogels. Prior to collecting data for statistic averaging, two interacting nanogels where generated and positioned as follows. First, one isolated nanogel was generated (and thermalized) applying the procedure described in a previous paper.45 Then a similar nanogel was generated replicating the first one and rotating by an angle with respect to it. Each simulation was repeated at six different initial relative orientations to minimize memory effects associated with such initial orientation. The second nanogel was originally positioned at a distance greater than their effective diameter to avoid overlapping and interactions. Figure 1a shows two nonoverlapping nanogels whose centers of mass (CM) are at a distance of 26 nm. Both nanogels (with 16 monomers per chains) are identical, but the second one has been rotated π/3 rad (around its CM). Finally, both nanogels were approached at a desired distance. Single-particle MC movements and translations of the whole nanogels were employed in this stage (approach) rejecting those movements that increased the center-to-center distance. Figure 1b shows the same nanogels after approaching. Their CMs are now at a distance of 20 nm. Simulation. After achieving the initial configuration (two nanogels whose centers of mass are at the desired distance), the system must be thermalized again. 1 × 107 MC steps were employed in each of the six thermalizations. After this process, data were collected for statistical averaging. 3 × 10 8 configurations were employed in the set of six simulations. During thermalization and averaging, small rotations of one of the nanogels were also employed (together with the singleubond(r ) =

MODEL AND SIMULATIONS

Model. The coarse-grained picture employed here is the socalled bead−spring model for polyelectrolyte, in which monomer units and ions are represented as spheres whereas the solvent is considered as a dielectric continuum, which has been extensively employed in the study of the collapse and adsorption of polyelectrolytes.28−37 In this survey, nanogels were simulated using a coarse-grained model quite similar to that employed by Claudio et al.,38 in which polymer chains are cross-linked in order to construct the polymer network. Four different ideal nanogels have been investigated with the same topology (illustrated in Figure 1). In all cases, the nanogel has 66 cross-linkers connected by 100 chains. The inner ones (26 out of 66) are connected by four polyelectrolyte chains. But the rest (the most external ones) are connected only by three or two chains. None of the nanogels have dangling chains. Thus, the main difference between the nanogels is the number of monomers per chain (Nm). Nanogels with Nm = 4, 8, 12, and 16 monomers per chain have been simulated here. Other properties, such as the total number of monomers per nanogel C


Article

Macromolecules particle MC movements and translations of the whole nanogel). Computation of the Effective Interaction Potential. For two nanogels at infinite dilution, the effective interaction potential, ueff(r), can be determined as ueff (r ) = −kBT ln(P(r )) + u0

(5)

where kB is Boltzmann’s constant, T is the absolute temperature, and P(r) is the normalized probability of finding their respective CMs at a separation r. The reader should also keep in mind that the interaction potential given by this equation is determined to within an additive constant (u0). In practice, however, this procedure is refined using the so-called umbrella sampling methods, which significantly improve the statistics of rarely visited configuration (for instance, strongly overlapping nanogels in our case) introducing a bias potential that favors such configurations. A particular and simple choice for our case is the following window potential: ⎧ 0 rmin < r < rmax w(r ) = ⎨ ⎩∞ otherwise

Figure 2. Effective interaction potential, ueff(r), obtained in the different windows from eq 5 for a nanogel with Nm = 4. Different colors are employed for each window.

the size of the nanogel through a length (R) that is proportional to the radius of gyration and can give us an idea of the geometrical radius of the nanogel considered in a first approximation as a sphere:

(6)

where rmin and rmax are the lower and upper limits of the window, respectively. Note that this infinite square well restricts the movements of the CM of both nanogels to distances between rmin and rmax, the window of r values which we want to focus on. The potential employed for the confinement in windows can have different functional forms (for instance, other authors use harmonic potentials). One of the main advantages of the infinite square well (eq 6) is that this potential does not require an explicit correction after sampling. In any case, some care must be taken when defining the positions and widths of the windows, particularly in the region where the effective potential is expected to vary steeply. Then, the width of the windows should be adjusted if sampling statistics is not good enough. In practice, a trial-and-error process is employed to determine the optimal width of this bias potential. Accordingly, the whole interval of r values is divided into a series of windows, applying eq 5 to compute ueff(r) from P(r) in each window. The complete functional form of ueff(r) can be reconstructed recalling that this interaction potential must be continuous in passing from one window to next. This condition can be fulfilled adjusting the additive constants of each window. In a previous work, this well-known procedure was employed to calculate the effective interaction energy between a negatively charged polyelectrolyte chain and a likecharged surface,46 which can be an instructive example of how such a procedure is applied for readers interested in further details. Figure 2 shows another example: The function ueff(r) obtained for a nanogel with Nm = 4. Different colors are employed for each window. In the window corresponding to the largest distances, u0 was adjusted so that ueff(r) tends to 0 at large r values. In the rest of windows, this additive constant was adjusted recalling that ueff(r) must be continuous.

R=

5/3 R g

(7)

This equation can be obtained computing the radius of gyration of a solid sphere of geometrical radius R (assuming that the mass is uniformly distributed) and solving for R. The roles play by R and Rg in the distribution of monomers (and mass) in the nanogel are discussed with the help of the Supporting Information. Here we only mention that the fraction of monomers included in a sphere of radius r = Rg is slightly greater than 0.50 and grows up to 0.9 for a distance r = R. In other words, R gives us a better idea of the geometrical dimension of our nanogels. Table 1 also displays the R values obtained from simulation for the systems studied in this work. The errors were computed from the fluctuations of this quantity during the runs corresponding to the generation of isolated nanogels.45 Interaction Potentials Obtained from Simulations. Figure 3 shows the effective interaction potential computed through eq 5 for these nanogels as a function of a normalized

■

RESULTS Size of Nanogels. As mentioned previously, one of the goals of this work is relating the nominal diameters appearing in expressions of ueff(r) with measurable quantities (computable quantities in this case). For polymer chains or networks, computer simulations can provide us the so-called radius of gyration, which can be easily obtained from static light scattering experiments. In this work, however, we will describe

Figure 3. Effective interaction potential, ueff(r), obtained for the different systems as a function of a center-to-center distance, r/2R. D


Article

Macromolecules center-to-center distance, r/2R. As can be seen, a considerable potential barrier appears between the nanogels in all the cases when their center-to-center distance is reduced. This behavior can be attributed to the high degree of cross-linking between polymer chains. Systems of highly penetrable coils, such as selfavoiding random walk polymers, were studied by Bolhuis et al.20 and Louis et al.47 These authors showed that the interactions between such polymers can be reasonably well represented by a Gaussian potential with a range of the order of the radius of gyration of individual coils. Apart from that, it should be stressed that each system exhibits a well-differentiated interaction potential. In other words, these functions cannot be described in terms of a master curve only depending on the normalized center-to-center distance. In this work, we will focus on the regime of slightly overlapping nanogels, where the differences between these curves are smaller. The data of Figure 3 have also been drawn in a log−log plot to see more clearly how the interaction potential decays at separations of the order of the particle diameter (Figure 4).

Table 2. Parameters of the Soft-Sphere Potential (Eq 1) for the System with Nm = 4, 8, 12, and 20 Describing Simulation Data in the Interval 0.95 < r/2R < 1.05a Nm 4 8 12 16

n 20.6 15.9 20.7 19.6

± ± ± ±

0.2 0.1 0.2 0.2

εn/kBT

correlation coefficient

± ± ± ±

−0.991 −0.994 −0.992 −0.989

0.295 0.398 0.1807 0.290

0.002 0.001 0.0007 0.002

a

These parameters were obtained plotting the logarithm of ueff(r)/kBT as function of the logarithm of r/2R. This figure is not explicitly shown in the paper because Figure 4 has the same aspect.

Figure 5. (ueff/kBT)0.4 against r for the different systems. Straight lines are linear fits for r/2R ≥ 0.85.

sketched points should be on a straight line whose slope and intercept are (εH/kBT)0.4/σH and (εH/kBT)0.4, respectively. As can be concluded from Figure 5, (ueff/kBT)0.4 does not exhibit a clear linear dependence in the whole interval of r values studied in this work. This is not surprising since the Hertzian potential is strictly valid for very small deformations. For that reason, linear fits were tried only in the regime of slight overlaps (more specifically, for r/2R ≥ 0.85). The slope (a), intersect (b), and correlation coefficient obtained by fitting this region through linear regressions are shown in Table 3. The correlation coefficients shown in this table suggest that the regime of slightly overlapping nanogels previously mentioned might be described by the Hertz potential. However, the comparison between the straight lines corresponding to the linear fits and simulation data reveals that the latter exhibit some slight curvature even fo5 r/2R ≥ 0.85, which suggests that this classical potential is just a first approach to the problem that should be improved. However, let us analyze the results of the linear regressions performed in depth. The parameters σH and εH can be computed from the slope and the intercept as σH = −b/a and εH/kBT = b2.5, respectively. These parameters are also included in Table 3. In order to establish a relationship between σH and 2R (the geometrical diameter of the nanogels seen as a sphere), it is quite instructive to plot the former versus the latter (Figure 6). Such a plot clearly shows a linear dependence between both parameters, which can be determined from a linear regression: σH ≅ (1.10 ± 0.05)2R. This obviously means that when two nanogels slightly overlap, the interaction can be described in terms of a diameter which is a bit greater than their geometrical diameter

Figure 4. Effective interaction potential, ueff(r), obtained for the different systems as a function of a center-to-center distance, r/2R, in a log−log plot. Straight lines represent the fits obtained from eq 1 in the interval of small deformations 0.95 < r/2R < 1.05.

This information can be useful for many applications in which nanogels come into contact by increasing their concentration but they only undergo small deformations. On the other hand, Figure 4 allows us to perform a preliminary analysis of the validity of the soft-sphere potential (eq 1). In a log−log plot, this functional form would look like a straight line of slope n, but data do not show this appearance. However, we could restrict ourselves to a narrower region. More specifically, the data in the region 0.95 < r/2R < 1.05 were fitted using eq 1 (assuming that σn = 2R). The parameters obtained from these fits are shown in Table 2. As can be seen, the n values describing the stiffness of our nanogels in the region of small deformations are of the order of 20. Paulin et al. obtained similar exponents for poly(methyl methacrylate) microgels,48 but Richtering and co-workers reported smaller values.7,8 Testing the Hertz Potential. As mentioned in the Introduction, the Hertz model has been widely applied, probably due to the lack of nonlinear treatments. Thus, it would be interesting to find out to what extent the Hertz potential (eq 2) can describe the effective potentials obtained in our simulations. To carry out this task, (ueff/kBT)0.4 was plotted against r (Figure 5). If the Hertz potential were valid, the E


Article

Macromolecules

Table 3. Parameters Obtained from the Linear Regressions of (ueff/kBT)0.4 against r and Physical Parameters Computed from Them no. of monomers per chain, Nm 4 8 12 16

a −0.461 −0.276 −0.243 −0.210

± ± ± ±

b 0.003 0.002 0.001 0.001

6.15 5.41 5.72 5.88

± ± ± ±

correlation coefficient −0.997 23 −0.995 43 −0.996 33 −0.997 97

0.03 0.03 0.03 0.03

σH (nm) 13.4 19.6 23.6 28.0

± ± ± ±

0.2 0.3 0.2 0.3

εH/kBT 94 68 78 84

± ± ± ±

9 7 8 8

Estimation of εH from Elastic Properties of Gels. As the Hertzian model has been extensively studied, one can even find theoretical expressions for εH. Indeed, this is a desirable feature in an effective potential when its parameters can considerably vary from one system to another. Thus, it would be quite instructive to derive theoretical estimates of εH following a different route and compare with our previous results. It is known that this parameter can be obtained as49 εH =

2YσH 3 15(1 − ν 2)

(8)

Here Y and ν are the so-called Young’s modulus and Poisson’s ratio, respectively. In turn, these two properties can be expressed as49,50

Figure 6. σH vs 2R. The straight line represents the linear fit.

Y = 3K (1 − 2ν)

(considered as spheres). This is somewhat logical since the nanogel does not have a well-defined and smooth surface of radius R, and consequently, its range (characterized by σH in the Hertzian potential) should be a bit larger than R. At this point, let us pay attention to the other parameter appearing in the Hertz potential, εH. Figure 7 displays this

ν=

3K − 2μ 2(3K + μ)

(9)

(10)

Here, K and μ are the bulk modulus and the shear modulus, respectively. Consequently, εH could be estimated from eq 8 if K and μ are previously known as follows. First, the values of these two parameters should be replaced in eq 10 to obtain Poisson’s ratio. Then Young’s modulus can be obtained from the values of K and ν through eq 9. Finally, eq 8 is applied using the values of Y and ν. In this work, K and μ were estimated from simulation of gels with the same number of monomers per chain. K can be computed from its definition as

⎛ ∂Π ⎞ K = φ⎜ ⎟ ⎝ ∂φ ⎠T

(11)

where φ is the polymer volume fraction of monomers in the polymeric gel and Π is the osmotic pressure. According to this expression, the osmotic pressure was obtained by means of coarse-grained simulations of gels following the procedure described elsewhere.51 The osmotic pressure was calculated for several φ values near the equilibrium polymer volume fraction (φe) of a gel with the same Nm value. Figure 8 shows the osmotic pressures as a function of φ (for φ values around φe) for the gels with Nm = 4, 8, 12, and 16. As can be seen, the osmotic pressure grows with the volume fraction of the polymer network, as expected, and this increase is more rapid for gels with short chains. In the case of polymer networks without external osmotic pressure, the equilibrium polymer volume fraction satisfies that Π(φe) = 0, which allows us to calculate φe. Equation 11 must be evaluated at φe in our case. For that reason the data in the neighborhood of φe were fitted by a linear function. The slope of the corresponding fits provides us the value of the derivative in eq 11 and the value of K. The other parameter required in eq 10 is the shear modulus, which was estimated from50

Figure 7. εH vs Nm. Squares stand for the results obtained from εH/ kBT = b2.5. Circles and triangles denote the estimates from Young’s modulus and Poisson’s ratio (see text for further details). Some data are slightly shifted in the x-axis so that the error bars can be clearly distinguished.

parameter (obtained previously from the fits shown in Figure 4) as a function of Nm. The error bars were estimated analyzing how εH varies when the range of r values in which the fits are performed is slightly modified, turning out to be of the order of 10% when the lower limit of this interval changes from 0.83·2R to 0.90·2R. Even admitting this uncertainty, Figure 7 suggests that the behavior of εH might not be monotonic. This possibility will be discussed later. F


Article

Macromolecules

Figure 8. Osmotic pressures as a function of φ for gels with Nm = 4, 8, 12, and 16. Solid lines are linear regressions in the neighborhood of φe.

⎛ φ ⎞1/3 1 μ = nckBT ⎜⎜ e ⎟⎟ 2 ⎝ φ0 ⎠

Figure 9. ⟨Ree2⟩/((Nm + 1)⟨Rbb2⟩) as a function of the polymer volume fraction. The arrows point to φ0.

adjustable parameter was employed in the estimates of εH from gels (ii) both procedures lead to values of similar order, which means that eqs 8−11 can provide us a first estimate of the εH values for nanogels. As mentioned earlier, nc was estimated from the actual number of chain in the simulation box at the reference polymer volume fraction, φ0. However, some authors claim that an effective nc value (rather than the real one) should be used.52 Others state that this parameter is related to the number of elastically active network chains (rather than all network chains).53 Thus, a different estimate of nc was carried out from the Flory−Rehner theory for gels. According to this formalism52

(12)

In this expression φ0 denotes a reference state in which the polymer chains describe random walks, and nc is concentration of chains just at such reference state. Both parameters were again estimated from simulations of gels with the same number of monomers per chain as follows. It is well established that the mean-square end-to-end distance ⟨Ree2⟩ of a random walk of N steps is given by ⟨R ee 2⟩ = N ⟨R bb 2⟩

(13)

where ⟨Rbb ⟩ is the mean-square bead-to-bead distance (measured from center to center). The reader should note that in evaluating ⟨Ree2⟩, the nodes are considered the ends of the chains in the network. Thus, the random walk associated with each chain contains Nm + 1 steps. Consequently 2

⟨R ee 2⟩ = (Nm + 1)⟨R bb 2⟩

⎡ ⎛ φ ⎞1/3⎤ ncvs ⎢ 1 φe φe + ln(1 − φe) + χφe = − ⎜⎜ e ⎟⎟ ⎥ NA ⎢⎣ 2 φ0 ⎝ φ0 ⎠ ⎥⎦ 2

(15)

In this expression, vs is the molar volume of the solvent and χ is the so-called polymer−solvent interaction (or Flory’s mixing) parameter. Given that solvent−solvent interactions, solvent− polymer interactions, or polymer−polymer attractive interactions (such as the hydrophobic forces or hydrogen bonds) are not explicitly considered in our coarse-grained simulations, χ = 0 was assumed, an effective nc value was derived from eq 15, and a new estimate for εH was obtained and also plotted in Figure 6 (estimate 2). As can be seen, these estimates are again of the order of the values obtained from the interaction potential, but the agreement is better. Equation 8 also reveals that εH depends on three parameters (Y, σH, and ν), whose behavior might exhibit opposite trends in some cases. For instance, σH grows with the number of monomers per bead. However, our calculations show that Young’s modulus decreases. Consequently, εH might exhibit a rich (and nonmonotonic) behavior with Nm (as mentioned above). Unfortunately, the large error bars of our estimates for εH do not allow us to quantitatively justify the existence of a minimum in εH. Modifying the Original Hertz Model. As mentioned in the Introduction, some authors have shown that the forces measured in some experiments could be fitted by the Hertz model if the exponent appearing in the corresponding potential and force expressions differs from the classical one.18 In force

(14)

The φ value satisfying eq 14 can be identified as φ0. Figure 9 shows the quotient ⟨Ree2⟩/((Nm + 1)⟨Rbb2⟩) as a function of the polymer volume fraction. This quotient is a decreasing function because the mean distance between cross-linkers (endto-end distance) reduces when the polymer volume fractions increases and the room for them becomes smaller. When this quotient is greater than 1, the polymer chains are larger than a random walk of the same number of steps, which is the reference for φ0. On the other hand, polymer chains are shorter than the random walk if this quotient is smaller than 1. The arrows point to φ0. Having determined this parameter, nc can also be estimated from simulations as nc ≅ 6φ0/πd3Nm (since 6φ0/πd3 is the number density of beads in the simulation cell). The εH values derived from eqs 8−12 are also plotted in Figure 7 for comparison (estimate 1). Their error bars were estimated assuming linear propagation of errors considering only the most important sources of uncertainty (the computation of the derivative in eq 11 and the determination of φe and σH). Because of the considerable number of mathematical operations involved, the results are affected by a significant error. At first sight, one might conclude that the values estimated from eq 8 are a bit greater than those obtained by simulations. In any case, it should be stressed that (i) none G


Article

Macromolecules

in the regime of small deformations the parameters characterizing the Hertzian potential can even be theoretically predicted to some extent: the effective diameter (its range) is a bit greater (10−15%) than the geometrical diameter whereas the interaction strength can be calculated from the effective diameter and some elastic properties (bulk modulus, shear modulus, Young’s modulus, Poisson’s ratio). In this work, such properties were computed from simulations of gels with the same number of monomers per chain. Our results also suggest that a modified Hertz potential can describe the interaction between the ideal nanogels studied in this work.

experiments, this exponent would be in the interval [1.5, 2.5], which is equivalent to [2.5, 3.5] in terms of potential energy. Inspired by this idea, the fitting procedure used for the original model was repeated now with ⎧ εH(1 − r /σH)m r < σH ueff (r ) = ⎨ r ≥ σH ⎩0 ⎪

⎪

(16)

where m ranges from 2 to 5. A maximum was found in the correlation coefficient of the linear regression performed for m values ranging from 3.5 to 4 (see Supporting Information), just at the upper limit of interval mentioned above. The figure for m = 3.5 analogous to Figure 6 (for for m = 2.5) show that linearity in the regime of small deformations studied in this work has improved a bit (see Supporting Information). In this case, the effective (or nominal diameter) is related to the geometrical one through σH ≅ (1.16 ± 0.06)2R. This suggests that the effective diameter is again a bit larger (16%) than 2R. In relation to this, we should keep mind that according to the structural data provided as Supporting Information, 99% of the monomers are inside a sphere whose diameter is 1.16·2R, in clear agreement with effective diameter obtained from the effective potential. Having determined the values of σH for each nanogel from the linear regressions mentioned in the previous paragraph, the data shown in Figure 3 can be represented as a function of 1 − r/σH in a log−log plot (Figure 10). As can be seen, the data

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.6b02333. Further information about the distribution of the monomers in the nanogel particles (local number density of monomers and number of monomers enclosed by a sphere of a given radius) and the fits obtained modifying the exponent of the Hertz model (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected] (M.Q.-P.). ORCID

Silvia Ahualli: 0000-0002-7329-0817 Manuel Quesada-Pérez: 0000-0003-0519-7845 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank the financial support from the following institutions: (i) “Ministerio de Economiá y Competitividad, ́ Plan Estatal de InvestigaciónCientifica y Técnica y de Innovación 2013-2016”, Projects FIS2016-80087-C2-1-P, FIS2016-80087-C2-2-P, and FIS2016-81924-REDT; (ii) European Regional Development Fund (ERDF).

■

Figure 10. Effective interaction potential obtained for the different systems as a function of 1 − r/σH in a log−log plot.

REFERENCES

(1) Ramos, J.; Forcada, J.; Hidalgo-Alvarez, R. Cationic Polymer Nanoparticles and Nanogels: From Synthesis to Biotechnological Applications. Chem. Rev. 2014, 114 (1), 367−428. (2) Bonham, J. A.; Faers, M. A.; van Duijneveldt, J. S. Non-Aqueous Microgel Particles: Synthesis, Properties and Applications. Soft Matter 2014, 10 (47), 9384−9398. (3) Heyes, D. M.; Brańka, a. C. Interactions between Microgel Particles. Soft Matter 2009, 5 (14), 2681. (4) Royall, C. P.; Poon, W. C. K.; Weeks, E. R. In Search of Colloidal Hard Spheres. Soft Matter 2013, 9 (1), 17−27. (5) Paloli, D.; Mohanty, P. S.; Crassous, J. J.; Zaccarelli, E.; Schurtenberger, P. Fluid-Solid Transitions in Soft-Repulsive Colloids. Soft Matter 2013, 9 (11), 3000−3004. (6) Mohanty, P. S.; Paloli, D.; Crassous, J. J.; Zaccarelli, E.; Schurtenberger, P. Effective Interactions between Soft-Repulsive Colloids: Experiments, Theory, and Simulations. J. Chem. Phys. 2014, 140 (9), 094901. (7) Senff, H.; Richtering, W. Temperature Sensitive Microgel Suspensions: Colloidal Phase Behavior and Rheology of Soft Spheres. J. Chem. Phys. 1999, 111 (4), 1705−1711. (8) Pyett, S.; Richtering, W. Structures and Dynamics of Thermosensitive Microgel Suspensions Studied with Three-Dimen-

corresponding to different systems collapse over a large region and exhibit a linear trend in a log−log graph. This suggests that the modified Hertz model given by eq 16 provides a good description of the potential of this sort of nanoparticles.

■

CONCLUSIONS In this work, the interaction potential between ideal neutral nanogels was computed through Monte Carlo coarse-grained simulations. A set of four nanogels with different numbers of monomers per chain was employed in such simulations. Our results reveal that neither the Hertzian potential nor the softsphere potential can describe the interaction between nanogels in the whole interval of the center-to-center distances studied. As a first approach, however, both might describe the interaction potential in the regime of slight overlap, but the description of the soft-sphere potential is a bit worse and restricted to a narrower region. Our simulations also show that H


Article

Macromolecules sional Cross-Correlated Light Scattering. J. Chem. Phys. 2005, 122 (3), 034709. (9) Wu, J.; Huang, G.; Hu, Z. Interparticle Potential and the Phase Behavior of Temperature-Sensitive Microgel Dispersions. Macromolecules 2003, 36 (2), 440−448. (10) Dintwa, E.; Tijskens, E.; Ramon, H. On the Accuracy of the Hertz Model to Describe the Normal Contact of Soft Elastic Spheres. Granular Matter 2008, 10 (3), 209−221. (11) Liu, K.-K. Deformation Behaviour of Soft Particles: A Review. J. Phys. D: Appl. Phys. 2006, 39 (11), R189−R199. (12) Brochu, H.; Vermette, P. Young’s Moduli of Surface-Bound Liposomes by Atomic Force Microscopy Force Measurements. Langmuir 2008, 24 (5), 2009−2014. (13) Hashmi, S. M.; Dufresne, E. R. Mechanical Properties of Individual Microgel Particles through the Deswelling Transition. Soft Matter 2009, 5, 3682−3688. (14) Schmidt, S.; Zeiser, M.; Hellweg, T.; Duschl, C.; Fery, A.; Möhwald, H. Adhesion and Mechanical Properties of PNIPAM Microgel Films and Their Potential Use as Switchable Cell Culture Substrates. Adv. Funct. Mater. 2010, 20 (19), 3235−3243. (15) Lin, D. C.; Shreiber, D. I.; Dimitriadis, E. K.; Horkay, F. Spherical Indentation of Soft Matter beyond the Hertzian Regime: Numerical and Experimental Validation of Hyperelastic Models. Biomech. Model. Mechanobiol. 2009, 8 (5), 345. (16) Seth, J. R.; Mohan, L.; Locatelli-Champagne, C.; Cloitre, M.; Bonnecaze, R. T. A Micromechanical Model to Predict the Flow of Soft Particle Glasses. Nat. Mater. 2011, 10 (11), 838−843. (17) Burmistrova, A.; Richter, M.; Uzum, C.; Klitzing, R. V. Effect of Cross-Linker Density of P(NIPAM-Co-AAc) Microgels at Solid Surfaces on the Swelling/shrinking Behaviour and the Young’s Modulus. Colloid Polym. Sci. 2011, 289 (5−6), 613−624. (18) Dimitriadis, E. K.; Horkay, F.; Maresca, J.; Kachar, B.; Chadwick, R. S. Determination of Elastic Moduli of Thin Layers of Soft Material Using the Atomic Force Microscope. Biophys. J. 2002, 82 (5), 2798−2810. (19) Holmqvist, P.; Mohanty, P. S.; Nägele, G.; Schurtenberger, P.; Heinen, M. Structure and Dynamics of Loosely Cross-Linked Ionic Microgel Dispersions in the Fluid Regime. Phys. Rev. Lett. 2012, 109 (4), 1−5. (20) Bolhuis, P. G.; Louis, A. A.; Hansen, J. P.; Meijer, E. J. Accurate Effective Pair Potentials for Polymer Solutions. J. Chem. Phys. 2001, 114 (9), 4296−4311. (21) Xu, X.; Kanduč, M.; Wu, J.; Dzubiella, J. Potential of Mean Force and Transient States in Polyelectrolyte Pair Complexation. J. Chem. Phys. 2016, 145 (3), 034901. (22) Narros, A.; Moreno, A. J.; Likos, C. N. Influence of Topology on Effective Potentials: Coarse-Graining Ring Polymers. Soft Matter 2010, 6 (11), 2435. (23) Narros, A.; Moreno, A. J.; Likos, C. N. Effects of Knots on Ring Polymers in Solvents of Varying Quality. Macromolecules 2013, 46 (9), 3654−3668. (24) Narros, A.; Moreno, A. J.; Likos, C. N. Architecture-Induced Size Asymmetry and Effective Interactions of Ring Polymers: Simulation and Theory. Macromolecules 2013, 46 (23), 9437−9445. (25) Truzzolillo, D.; Bordi, F.; Sciortino, F.; Sennato, S. Interaction between like-Charged Polyelectrolyte-Colloid Complexes in Electrolyte Solutions: A Monte Carlo Simulation Study in the Debye-H??ckel Approximation. J. Chem. Phys. 2010, 133 (2), 024901. (26) Verso, F. Lo; Yelash, L.; Egorov, S. A.; Binder, K. Interactions between Polymer Brush-Coated Spherical Nanoparticles: The Good Solvent Case. J. Chem. Phys. 2011, 135 (21), 214902. (27) Marzi, D.; Likos, C. N.; Capone, B. Coarse Graining of StarPolymer - Colloid Nanocomposites. J. Chem. Phys. 2012, 137 (1), 014902. (28) Ulrich, S.; Seijo, M.; Stoll, S. The Many Facets of Polyelectrolytes and Oppositely Charged Macroions Complex Formation. Curr. Opin. Colloid Interface Sci. 2006, 11 (5), 268−272.

(29) Dobrynin, A. V. Theory and Simulations of Charged Polymers: From Solution Properties to Polymeric Nanomaterials. Curr. Opin. Colloid Interface Sci. 2008, 13 (6), 376−388. (30) Dias, R. S.; Pais, A. Polyelectrolyte Condensation in Bulk, at Surfaces, and under Confinement. Adv. Colloid Interface Sci. 2010, 158 (1−2), 48−62. (31) Jorge, A. F.; Sarraguca, J. M. G.; Dias, R. S.; Pais, A. Polyelectrolyte Compaction by pH-Responsive Agents. Phys. Chem. Chem. Phys. 2009, 11 (46), 10890−10898. (32) de Carvalho, S. J.; Metzler, R.; Cherstvy, A. G. Critical Adsorption of Polyelectrolytes onto Charged Janus Nanospheres. Phys. Chem. Chem. Phys. 2014, 16 (29), 15539−15550. (33) de Carvalho, S. J.; Metzler, R.; Cherstvy, A. G. Inverted Critical Adsorption of Polyelectrolytes in Confinement. Soft Matter 2015, 11 (22), 4430−4443. (34) Clavier, A.; Seijo, M.; Carnal, F.; Stoll, S. Surface Charging Behavior of Nanoparticles by Considering Site Distribution and Density, Dielectric Constant and pH Changes - a Monte Carlo Approach. Phys. Chem. Chem. Phys. 2015, 17 (6), 4346−4353. (35) Stoll, S. In Soft Nanoparticles for Biomedical Applications; The Royal Society of Chemistry: 2014; Chapter 10, pp 342−371. (36) Carnal, F.; Clavier, A.; Stoll, S. Modelling the Interaction Processes between Nanoparticles and Biomacromolecules of Variable Hydrophobicity: Monte Carlo Simulations. Environ. Sci.: Nano 2015, 2 (4), 327−339. (37) Uhlik, F.; Kosovan, P.; Zhulina, E. B.; Borisov, O. V. ChargeControlled Nano-Structuring in Partially Collapsed Star-Shaped Macromolecules. Soft Matter 2016, 12 (21), 4846−4852. (38) Claudio, G. C.; Kremer, K.; Holm, C. Comparison of a Hydrogel Model to the Poisson-Boltzmann Cell Model. J. Chem. Phys. 2009, 131 (9), 094903. (39) Mann, B. A.; Everaers, R.; Holm, C.; Kremer, K. Scaling in Polyelectrolyte Networks. Europhys. Lett. 2004, 67 (5), 786−792. (40) Mann, B. A.; Holm, C.; Kremer, K. Swelling of Polyelectrolyte Networks. J. Chem. Phys. 2005, 122 (15), 154903. (41) Yin, D. W.; Yan, Q. L.; de Pablo, J. J. Molecular Dynamics Simulation of Discontinuous Volume Phase Transitions in HighlyCharged Crosslinked Polyelectrolyte Networks with Explicit Counterions in Good Solvent. J. Chem. Phys. 2005, 123 (17), 174909. (42) Yin, D.-W.; Olvera de la Cruz, M.; de Pablo, J. J. Swelling and Collapse of Polyelectrolyte Gels in Equilibrium with Monovalent and Divalent Electrolyte Solutions. J. Chem. Phys. 2009, 131 (19), 194907. (43) Kosovan, P.; Richter, T.; Holm, C. Modeling of Polyelectrolyte Gels in Equilibrium with Salt Solutions. Macromolecules 2015, 48 (20), 7698−7708. (44) Ghosh, S. K.; Cherstvy, A. G.; Metzler, R. Deformation Propagation in Responsive Polymer Network Films. J. Chem. Phys. 2014, 141 (7), 074903. (45) Quesada-Perez, M.; Martin-Molina, A. Monte Carlo Simulation of Thermo-Responsive Charged Nanogels in Salt-Free Solutions. Soft Matter 2013, 9 (29), 7086−7094. (46) Luque-Caballero, G.; Martín-Molina, A.; Quesada-Pérez, M. Polyelectrolyte Adsorption onto like-Charged Surfaces Mediated by Trivalent Counterions: A Monte Carlo Simulation Study. J. Chem. Phys. 2014, 140 (17), 174701. (47) Louis, A. A.; Bolhuis, P. G.; Hansen, J. P. Mean-Field Fluid Behavior of the Gaussian Core Model. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2000, 62 (6), 7961−7972. (48) Paulin, S. E.; Ackerson, B. J.; Wolfe, M. S. Equilibrium and Shear Induced Nonequilibrium Phase Behavior of {PMMA} Microgel Spheres. J. Colloid Interface Sci. 1996, 178 (1), 251−262. (49) Riest, B. J.; Mohanty, P.; Schurtenberger, P.; Likos, C. N. Coarse-Graining of Ionic Microgels: Theory and Experiment. Z. Phys. Chem. 2012, 226, 711−735. (50) Hirotsu, S. Softening of Bulk Modulus and Negative Poisson’s Ratio near the Volume Phase Transition of Polymer Gels. J. Chem. Phys. 1991, 94 (5), 3949−3957. (51) Quesada-Pérez, M.; Ahualli, S.; Martín-Molina, A. ThermoResponsive Gels in the Presence of Monovalent Salt at Physiological I


Article

Macromolecules Concentrations: A Monte Carlo Simulation Study. J. Polym. Sci., Part B: Polym. Phys. 2014, 52 (21), 1403−1411. (52) Quesada-Perez, M.; Alberto Maroto-Centeno, J.; Forcada, J.; Hidalgo-Alvarez, R. Gel Swelling Theories: The Classical Formalism and Recent Approaches. Soft Matter 2011, 7 (22), 10536−10547. (53) Dusek, K.; Duskova-Smrckova, M. Network Structure Formation during Crosslinking of Organic Coating Systems. Prog. Polym. Sci. 2000, 25, 1215−1260.

J


Interaction between Ideal Neutral Nanogels: A Monte Carlo Simulation

Recommend Documents