A New Self-Consistent Field Model of Polymer/Nanoparticle Mixture

Field-theoretical method is efficient in predicting assembling structures of polymeric systems. However, it’s challenging to generalize this method to study the polymer/nanoparticle mixture due to its multi-scale nature. Here, we develop a new field-based model which unifies the nanoparticle description with the polymer field within the self-consistent field theory. Instead of being “ensemble-averaged” continuous distribution, the particle density in the final morphology can represent individual particles located at preferred positions. The discreteness of particle density allows our model to properly address the polymer-particle interface and the excluded-volume interaction. We use this model to study the simplest system of nanoparticles immersed in the dense homopolymer solution. The flexibility of tuning the interfacial details allows our model to capture the rich phenomena such as bridging aggregation and depletion attraction. Insights are obtained on the enthalpic and/or entropic origin of the structural variation due to the competition between depletion and interfacial interaction. This approach is readily extendable to the study of more complex polymer-based nanocomposites or biology-related systems, such as dendrimer/drug encapsulation and membrane/particle assembly.

Scientific RepoRts | 6:20355 | DOI: 10.1038/srep20355 particle density in our model is eventually discrete and can represent individual particles as in the cavity model and other particle-based methods. This kind of particle density is no longer an "ensemble-averaged" distribution but of instantaneous nature. The morphology we predict reflects the ensemble-averaged density distribution of polymers with individual nanoparticles located at the "most probable" positions. The motivation is three-fold: 1) patterns with instantaneous particle locations can better exhibit delicate structures (especially when the particles are packed) and be compared with experimental results; 2) recovering the particle coordinates allows proper treatment on the EV interactions and the polymer-particle interface; 3) due to 2), the evaluation of the polymer-particle interfacial energy is more appropriate and the contribution of the depletion effect can also be involved. We use this model to study the simplest system of nanoparticles immersed in dense homopolymer solution. We focus on the structural variations under diverse polymer-particle interfacial interactions and depletion thickness.

Model and Method
We consider a system consisting of a mixture of n p homopolymers and n c spherical solid particles. We use the grand canonical ensemble to describe the polymers, i.e. n p is determined by the chemical potential and may vary under different situations. While n c is assigned and fixed in the calculation. The radius of all particles is R c . The coarse-grained polymer chain is composed of N segments of size σ. The bonded interactions of Gaussian chain are quantified by an elastic potential energy: represents the configuration of the polymer chain, R g 0 is the unperturbed radius of gyration. We employ the quadratic compressible model to address the EV interaction between segments in the presence of particles: is the number density of particles. The surface profile of particles is smoothed by the hyperbolic tangent function for numerical efficiency. Note that, for conciseness, we directly use the "ensemble-averaged" symbols in the above and following energy expressions, since they are interchangeable with instantaneous ones in the mean field approximation.
In the calculation, the particle density gradually evolves into individual particles. And in this process it's important to decide whether a location is in the hardcore or non-hardcore regions of particles. Hence, we introduce φ hc in Eq. (2) as a threshold quantity. We set φ φ = .
( ) max r 0 3 [ ] hc r c   , i.e., about one third of the maximum particle concentration in the system. In the hardcore (non-hardcore) region, ). φ hc varies during the calculation and eventually, φ = . (2) only serve as a judgment to pick the non-hardcore regions or the hardcore regions where the total concentration is less than one. In these two cases, only the weak polymer-polymer EV interaction characterized by κ −1 is considered. Otherwise, the strong polymer-particle EV interaction characterized by κ − h 1 (Eq. (6)) is triggered to ensure the total concentration in the hardcore region does not exceed 1.
(2) is introduced to take into account the depletion layer surrounding the surface of each particle and the capability of the layer to overlap with each other and with the entity of particles. This capability of overlapping leads to the so-called depletion effect (see Fig. 1). We approximate that the density profile of polymers in the depletion layer is described by a hyperbolic tangent function in the case of neutral polymer-particle interfacial interaction. Equivalently, we set an effective particle concentration (not real) to mimic the depletion layer felt by polymers. The particle at position ′ r  generates an effective particle concentration at position r  : where l is the grid size of the simulation box. The thickness of depletion layer, ξ D , is related to the rigidity of the polymer chain 43,44 . The Heaviside function in Eq. (4) guarantees the position r  is outside the entity of the particle at position  ′ r . The effective particle concentration is then: The first max function addresses the overlap between the depletion layer and the entity of particles, while the second one addresses the overlap between the depletion layers of nearby particles. Figure 1 shows the concentration plots of two contact complete particles in our model. Apparently, the particles are surrounded by a depletion layer mimicked by effective particle concentration felt by polymers. And, the depletion layers can overlap (not summation) with each other and with the entity of particles. The strong polymer-particle EV interaction is triggered where the total concentration in the hardcore region of particles is larger than one (the first term of Eq. (6)): The second term of Eq. (6) represents particle-particle hardcore repulsion which is set to be proportional to the overlap "volume", V op , between two nearby particles. κ − h 1 is chosen to be large enough to avoid overlap between particles and between polymer and particle. This expression for particle-particle EV repulsion is appropriate only in terms of instantaneous density, which is consistent with the instantaneous nature of the particle density in our model. Note that the ensemble-averaged particle concentrations do overlap.
The chemical nature of the polymer and particle is encoded in an interfacial energy of exponential form 25 : χ is the dimensionless strength of the interfacial interaction or relative affinity between polymer and particle (positive for repulsion and negative for attraction). ∆ denotes the spatial range. Besides the above potentials of real origins, an artificial double-well like potential is introduced to force the formation of individual particles: λ is the strength of the two harmonic potentials which drive the particle number density at each grid in the simulation box toward 0 or − l d . This artificial potential embodies the inseparability of a real particle. Note that this potential is position-unbiased, i.e., it does not directly influence the spatial arrangements of particles.
We fix σ = . R 4 08 g 0 . The chemical potential of polymer, µ, is chosen so that the concentration of the bulk polymer is 1. We set κ = . − 3 33 1 , which corresponds to the compressibility of polymethylmethacrylate melt at 450 K 45,46 . The SCFT calculations are performed in two dimensions (see Supplementary Information for the full SCFT equations) 47,48 . As in the real-space screening method 48 , the calculation is started with randomly generated potential fields for both the polymer and particle. Annealing process (gradually increasing λ) is used for the gradual formation of the individual particles and for effectively searching the morphology with low energy. Rather than predicting the configurations corresponding to the global minima, we only look for the specific morphologies showing the features of particle assembly at certain conditions.

Results and Discussion
Assembly of small nanoparticles. We first consider the case of small nanoparticles ( The depletion effect is neglected (by setting ξ = 0 D ), since the size of the nanoparticle is comparable to that of the segment. The typical assembly morphologies of varying the strength of interfacial interaction are shown in Fig. 2, (a-j). The particle concentration in these final morphologies can represent individual particles showing clear interfaces with surrounding polymers. The packing or dispersion of these particles is directly visible. For weak repulsion, χ = .
0 5 (Fig. 2(a,b)), the particles macroscopically separate from the polymer matrix and aggregate contactly (closely packed) into one large cluster (Contact Aggregation, CA) to minimize the interfacial energy. It is a typical enthalpic-driven phase separation of a binary mixture. Non-uniform clusters of contact-aggregated particles (Contact-Aggregation Clusters, CAC) are formed for very small repulsion (or χ = 0, where weak entropic attraction is caused by the hyperbolic tangent surface profile of particles) (Fig. 2(c,d)). As expected 26 , in the case of weak attraction (Fig. 2(e,f)), particles are dispersed randomly in the polymer matrix, each one of which is coated with a layer of slightly more concentrated polymers (Random Dispersion, RD). As the attraction is enhanced, the bound polymer layers get denser and particles get closer to share the bound layers. This aggregation of particles is ascribed to the bridging effect 26 . Irregular domains rich in both polymers and particles (Loosely Bridging Aggregation, LBA) are formed ( Fig. 2(g,h)) at intermediate attraction strength. For strong attraction (χ = − . 3 0), the nanoparticles aggregate closely (Closely Bridging Aggregation, CBA). The system can be viewed as composed of two separated phases: bulk polymer solution and the concentrated phase rich in both polymers and particles ( Fig. 2(i,j)). Note that, intuitively, the energy of the intermediate morphologies ( Fig. (c,d), (g,h)) should not be the global minima. But we still consider them as featured morphologies and ascribe their emergence to the configurational entropy of particles. This entropy takes effect during the search of the energy landscape but is missed in the energy expression. Fig. 2, we introduce two useful quantities: the number of pairs of bridging-connected particles, B n and the mean distance between a particle and its six nearest neighbors, D 6 . B n reflects the degree of bridging aggregation of particles from the view of amount. D 6 reflects the average degree of packing of particles from the view of distance. Two particles are deemed bridging-connected if a) their surface-to-surface distance is less than σ 4 , and b) the concentration of polymer in between exceeds 1.2. The choice of σ 4 is based on the calculation of two-particle potential of mean force in the dilute particle limit, which shows the range of attraction well for χ < 0 extends to σ 4 . The number 1.2 (> 1) is empirical, which does not influence the trend of B n . The variation of these two quantities with χ is shown in Fig. 3(a). Every data point is averaged over 10 independent runs. For different initial seeds, the overall feature of the final morphologies is the same but the packing details of particles are different, which causes a weak fluctuation in the energy and entropy (see the error bars in Fig. 3(b)). When χ decreases from 0.5 to 0, D 6 grows rapidly due to the increase of the number of particles at the interface (corresponding to the variation of the morphology from CA to CAC). A sharp transition happens between χ = 0 and χ = − . 0 2 that particles are no longer in contact (i.e., the transition of the morphology from CAC to RD). For weak attraction ( χ − < < 1 0 ),  B 0 n and D 6 reaches a maximum plateau, which implies the well dispersion of particles in the polymer matrix. As the attraction is enhanced, B n (D 6 ) increases (decreases) gradually, indicating the formation of more and more bridging connections or the closer and closer aggregation of particles (i.e., the formation of LBA). At very large attraction, B n (D 6 ) saturates at ~300 pairs of bridging connections (average center-to-center distance of ~σ 4 ), corresponding to the formation of the CBA morphology. The calculation of B n and D 6 helps determine the morphology diagram in Fig. 3(c).

Morphology diagram and thermodynamic analysis. To quantitatively analyze the morphologies in
To thermodynamically understand the variation of the morphologies, we take the particle distribution in RD morphology, specifically the distribution of particles in Fig. 2(e), as the reference. We calculate the differences of grand potential (∆Ω), potential energy (∆E) and entropy per polymer chain (∆S s ) between the "equilibrium" and the reference morphology at various χ (Fig. 3(b) and see Supplementary Information for the expressions of Ω, E and S s ). For χ > 0, ∆E is negative and ∆S s is positive, i.e., both the enthalpy and entropy are the driving force for the formation of contact aggregation of particles. The values of ∆E and ∆Ω are close, which means enthalpy is the main driving force. For χ = 0, the potential energy does not change (∆ = E 0) with the rearrangement of particles, and the negativity of ∆Ω comes from the increase of entropy of polymer chains, i.e., the CAC morphology at Depletion vs. interfacial attraction in the case of large nanoparticle. We also investigate large nanoparticles immersed in the polymer matrix. The radius of particles is set to be σ = R 5 c (≈ . , R 1 22 g 0 ). Since the particle is much larger than the segment, we turn on the depletion effect, by controlling the parameter ξ D . Larger ξ D implies more rigid polymer chain 43,44 . Here, we focus on the competition between the depletion and the interfacial interaction. When χ = − . 0 2, the system is close to the boundary between RD and CA phases. Particle distributions for σ ∆ = .
1 5 D are shown in Fig. 4  (a-c), respectively. The particle distribution changes from RD to CA along the path of increasing ξ D at constant σ ∆ = .
0 5 . While, on the contrary, the distribution changes from CA to RD along the path of increasing ∆ at constant ξ σ = .
1 5 D . Hence, we have the conclusion that narrow depletion layer (or flexible polymer chains) and/ or "long"-range weak polymer-particle attraction can facilitate the dispersion of large nanoparticles into polymer matrix.
We take the particle distributions of Fig. 4(a,b) as the "standard" RD and CA distributions, respectively and calculate the differences of ∆Ω = Ω − Ω 0 5 and varying ∆ with constant ξ σ = .
1 5 D (Fig. 4(d,e)). Both ∆E and ∆S s at all cases are positive, i.e., the aggregation of particles is energetically unfavored but entropically favored (depletion effect). In Fig. 4(d), the transition from RD to CA happens between ξ σ = D and ξ σ = .
1 5 D where ∆Ω becomes zero. ∆S s (∆E) increases (decreases) with ξ D , implying that the increase of depletion thickness (rigidity of poly- mer chain) enhances the aggregation of particles not only through the well-known depletion entropic effect but also through the alleviation of the penalty of interfacial energy. The transition from CA to RD in Fig. 4(e) happens at ∆ slightly larger than σ . 0 5 . The decrease of ∆S s and the increase of ∆E imply that the contribution of enthalpy (entropy) in the structural variation is enhanced (weakened) with the increase of ∆.

Conclusion
We have introduced a field model of PNCs which realizes the discrete description of particles in the final morphology. It is able to predict the collective assembly of particles in the polymer matrix with careful consideration of the EV interactions, depletion effect and interfacial interaction. This model allows investigations of the bridging and depletion effects on the multi-particle level instead of calculations based on two-particle correlations 25,26 . It is also able to reveal the entropic and enthalpic contributions in the variation of morphologies. Overall, it is a valuable approach to explore and analyze the rich mesostructures or particle-induced defects 49 in polymer-based nanocomposites and is also readily extendable to the study of biology-related systems, such as dendrimer/drug encapsulation 50 and membrane/particle assembly 51 .