Effect of size and charge asymmetry on aggregation kinetics of oppositely charged nanoparticles

We report a theoretical and experimental study of the aggregation kinetics of oppositely charged nanoparticles. Kinetic Monte Carlo simulations are performed for symmetric, charge-asymmetric and size-asymmetric systems of oppositely charged nanoparticles. Simulation results show that both the weight and number average aggregate size kinetics exhibit power law scaling with different exponents for small and intermediate time of evolution. The qualitative behavior of the symmetric and the size asymmetric system are the same, but the charge asymmetric system shows anomalous behavior for intermediate to high particle concentrations. We also observe a strong dependence of power law exponents on the particle concentration. Radius of gyration of the cluster that indicates how nanoparticles inside a cluster are distributed around the center of mass of the cluster shows a non-monotonic time evolution with pronounced peak at higher particle concentration. The dependence of particle concentration on aggregation kinetics as observed by predictive numerical simulation is further verified experimentally by monitoring the time evolution of aggregate size of nanoparticles assemblies of Poly (methacrylic acid) (PMMA) nanoparticles functionalized with oppositely charged ligands. These size and charge tunable asymmetric polymeric nanoparticles were synthesized by modified miniemulsion technique. The integrated approach for studying nanoparticles aggregation as described here renders new insights into super structure formation and morphology optimization which can be potentially useful in the design of new materials, such as organic photovoltaics.

(2019) 9:3762 | https://doi.org/10.1038/s41598-019-40379-y www.nature.com/scientificreports www.nature.com/scientificreports/ predictive modeling. Several theoretical models 22,23,29,33 to explain aggregation of sub-micrometer spherical particles have been developed and experimentally validated, the most common being DLVO theory. Computer simulation plays a crucial role to understand the internal dynamics of the system and provide a complimentary tool for the rational design of chemical reactors for high-throughput synthesis of nanoparticles.
Monte Carlo(MC) simulation is very useful tool to study the final equilibrium structure of super lattice formed by aggregation of NPs and the effect of various interactions on these structures 12 . Using importance sampling of configuration space, MC simulation efficiently reaches the final equilibrium state by skipping several metastable states. But during the evolution to reach equilibrium, system passes through various non-equilibrium complex structures which are very different from equilibrium structures 34 . For instance, various studies have shown that system of oppositely charged nanoparticles are often trapped in non-equilibrium "gel"-like (percolated) structures at high packing fractions, and do not evolve to minimum energy crystalline structures due to slow system dynamics [35][36][37] . In one of our previous study, we have found that these kinetically trapped percolated structures formed by self-assembly of electron conducting and hole conducting polymer NPs are the bulk heterojunction morphologies for efficient organic solar cells 37 .
In this paper, using kinetic Monte Carlo (kMC) simulations, we study the effect of charge and size asymmetry on aggregation kinetics of oppositely charged NPs. kMC scheme used in this simulation has been used in various other studies [36][37][38][39] and is argued to be more efficient than Brownian dynamics simulation because it employs a longer time step without compromising the numerical stability. The predictive modeling is validated experimentally by synthesizing size and charge tunable functionalized polymeric NPs with modified mini-emulsion method. Polymeric NPs functionalized with oppositely charged surfactants are chosen as model system to study aggregation behavior as they can be easily synthesized with different sizes and charges by tuning polymer's molecular weight, polymer to surfactant ratio, sonication time, solvent evaporation rate etc. The time evolution of the aggregate size of these PMMA NPs is studied by DLS.

Methodology
Numerical Simulation (Model and Method). Implicit-solvent simulations are performed on the solution containing two types, type A (+) and type B (−), of oppositely charged spherical NPs of diameter σ A and σ B and valence z A and z B , respectively. The system containing total N T NPs satisfy the electroneutrality condition n A z A + n B z B = 0 and n A + n B = N T , where n A and n B are number of NPs of type A and type B respectively. NPs interact via electrostatic interactions, which are modeled using a pairwise screened Coulomb potential between i th and j th NP at distance r given by 40 The soft-core repulsion between NP is modeled using Lennard-Jones (LJ) potential given by 12 6 where i, j = A, B and σ ij = (σ i + σ j )/2 is the average diameter of i th and j th NPs. We define dimensionless Bjerrum length as λ πεε σ = e kT /4 B B 2 0 0 and Debye screening length as κ ε, ε 0 , k B , T and c s are the dielectric constant of solvent, permittivity in free space, Boltzmann constant, absolute temperature and salt concentration respectively. In simulation we have used truncated and shifted form of both electrostatic and soft-core potential to improve computational efficiency, that is, are cutoff distances for Columbic interaction and softcore LJ repulsion respectively. We have not included the Van der Walls interactions in our simulations because we are interested in studying the effect of electrostatic interactions on aggregation kinetics and strength of electrostatic interactions we are using are of the order of 80k B T.
Simulation algorithm comprises of following steps 38 : a) Random distribution of NPs inside cubic box of length , where η is the packing fraction, defined as the fraction of volume occupied by NPs. b) Trial displacement of the center of a randomly picked NP to the surface of a sphere of radius a around the center of NP. We choose a random point (θ cos (2 1) 1 ) on the surface of sphere where u and v are uniformly distributed random variables between (0, 1); θ and ϕ are the polar angle and azimuthal angle, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ c) Move is accepted or rejected using the Glauber transition probability = where ΔE is the energy change for the trial displacement.
One kMC sweep consist of N T trial moves. Time step Δt of a sweep is related to the step size a as Δ = t a D 12 2 , where D is the diffusion coefficient of NP. For NP of diameter σ 0 , characteristic diffusion time scale is . In theory, Δt is calculated by assuming the energy change in each sweep is small. In previous kMC study 39 , a = 0.02σ 0 is observed as the proper step size and we have used this value of a in our simulation. Note that the current approach do not explicitly account for changes in diffusion coefficient on aggregation, as elaborated in previous studies 36,38,39 . However, this limitation is not expected to be of much consequence in the current work, since we focus on early stage kinetics of NP aggregation when the aggregate size is relatively small except for the high particle concentration case where percolation occurs.

Experiment (Materials and Method).
Chloroform, SDS, CTAB, PMAA and acetone were purchased from Sigma-Aldrich. All the materials were used as obtained. Millipore water was used in all parts of experiment including NP synthesis and DLS study. We synthesized both positively and negatively charged PMAA NPs of different concentrations with mini-emulsion technique using anionic surfactant, SDS (sodium dodecyl sulfate) and cationic surfactant, CTAB (cetyltrimethyl ammonium bromide), respectively. In a typical fabrication process, SDS/CTAB and PMAA are separately dissolved in distilled water and chloroform, respectively. Polymer solution was heated to 50 °C to ensure complete solubility. The NPs were then formed by adding the polymer solution to the aqueous SDS/CTAB solution under probe sonication with constant heat environment at 50 °C and stirring at 1000 rpm by hot plate. After complete cycle of probe sonication for 6 minute, dispersion was kept on heating at 65 °C for 45 minute to remove chloroform. To perform aggregation kinetics experiment, we prepared two sets by mixing both positively and negatively charged NPs together. In set 1, 200 ml was prepared by mixing distilled water with 25 μl of each positive and negative NPs. Similarly in set 2, 50 μl of each positive and negative NPs were mixed in distilled water to get 200 ml solution. Average aggregate size(M N ) and zeta potential(ζ) of NPs was recorded by Malvern Nanosizer DLS instrument.

Results and Discussion
Simulation Results. We numerically studied the aggregation kinetics for the following three sets of oppositely charged NPs: 1) Symmetric NPs: Both type A and type B have same size (σ A = σ B ) and equal and opposite valence (z A = −z B ). 2) Charge Asymmetric NPs: Both type A and type B have same size (σ A = σ B ) but unequal valence (z A = −2.0z B ). 3) Size Asymmetric NPs: type A and type B have unequal size (σ A = 1.5σ B ) but equal and opposite valence (z A = −z B ).
Starting from random distribution of NPs inside the simulation box we computed the aggregate size distribution of NPs at different times for three different packing fractions (η = 0.005,0.025,0.06). Aggregate size of a cluster is calculated by computing center to center distance(r) between two NPs and two NPs are considered to be part of the same cluster if r < 1.5σ 0 . We computed mean average aggregate size ( = ∑ ∑ M n j n / N j j j j ), weight average aggregate size ( = ∑ M n j n j / w j j j 2 ) and mean radius of gyration square as a function of kMC steps, where the number of aggregates of size j is n j , N c (= ∑ n j j ) is total number of aggregates in the system.
is the radius of gyration square of i-th aggregate of size N i and center of mass We have chosen simulation parameters, λ B = 81, κ −1 = σ 0 such that we are working in the NPs aggregation regime 36 . Figure 1 shows the time evolution of M N and M w of all the three sets for three different packing fractions (η = 0.005 (), 0.025(•), 0.06(Δ)). We observed two different aggregation kinetics regimes and time dependence of average aggregate size can be written as for all three sets. The β 1 and β 2 values characterize aggregation rate at short and intermediate time respectively. Both β 1 and β 2 strongly depend on the concentration of the NPs in the system as evident from Fig. 1. For η = 0.06, β 1 (short time aggregation regime) was not captured due to initial fast clustering and only β 2 values are given in Fig. 1. However, initial kinetics regime is slightly visible for η = 0.06 in size asymmetric case where aggregation is slow due to large size of the NPs. As η increases, both β 1 and β 2 increases for all the three cases. But, there is no substantial difference in β 2 values for η = 0.025 and η = 0.005 in charge asymmetric case. Furthermore, for all the three sets, we observed β 2 > β 1 for all η, but again in charge asymmetric case this difference is not very predominant as η increases. The behavior of exponents β 1 and β 2 in charge asymmetric case for intermediate to high η hints towards anomalous behavior of charge asymmetric system. We further plotted M w vs t/τ 0 and observed similar time dependence behavior as M N vs t/τ 0 , therefore we write www.nature.com/scientificreports www.nature.com/scientificreports/ The increase in the exponent β 1 and β ′ 1 as η increases is due to decrease in the average separation between NPs. This results in an increase in the Columbic force between NPs and hence faster aggregation at initial time for large η values is observed. Strong Columbic attraction due to high charge on NPs leads to much faster initial aggregation in charge asymmetric case as compared to symmetric and size asymmetric case (see Fig. 2A,B, blue data points are above red for short time). The aggregation become even faster at intermediate time for all the cases except for larger η values in charge asymmetric case (see Fig. 1). At intermediate time there are two processes that will determine the growth of the aggregate -one is the strength of electrostatic interactions between the NP clusters present in the system that determine how fast the clusters agglomerate and another is the size of the clusters that are agglomerating. www.nature.com/scientificreports www.nature.com/scientificreports/ The interplay between these two processes will determine the β 2 and β ′ 2 behavior. At initial time, small aggregates mainly comprising of dimers or trimer starts forming (see Fig. 2C, t = 0.02), rate of which depend on the concentration of particles. In charge symmetric systems dimers are more stable whereas in charge asymmetric case trimers are more favorable due to neutralization of charge on the cluster. These dimers and trimer formed in the system will now interact with other small clusters present in the system which are mostly dimers and trimers. Dimers will behave like electric dipoles and trimer will behave like electric quadrupole 41 . The interaction between two dimers or trimers will depend on the electric field generated by these small clusters. For the sake of simplicity for a diluted solution we can assume that electric field ∝ E r 1 3 for dimers and ∝ E r 1 4 for trimers, where r is the distance of the point from the cluster where electric field is calculated. This shows that the interaction between dimers will be stronger than trimers and hence the agglomeration become slower in the case of asymmetric charge at intermediate time when compared to the symmetric case. We should emphasize here that the actual interaction will also depend on the orientation of clusters with respect to each other. Moreover, if the solution is not sufficiently dilute, other higher order terms like , r r 1 1 5 6 etc. cannot be neglected from the electric field expression. In Fig. 2B  (t = 0.02 snapshot), one NP of high valence (charge asymmetric case) is surrounded by many NPs of opposite charge which leads to more screening and excluded volume effect 42 at intermediate time. This leads to large reduction in the strength of electrostatic interaction between the clusters and hence slower agglomeration of clusters. The increase in the size of agglomerated cluster due to the merging of smaller clusters is dominated by slower agglomeration of clusters in charge asymmetric case despite the fact that the asymmetric case have bigger clusters to agglomerate as compared to the symmteric case. Thus, we observe a crossing of aggregate size between symmteric and charge asymmetric case (see Fig. 2A,B). We further calculated polydispersity index (PDI) which is defined as the ratio of weight average aggregate size and number average aggregate size and given by PDI = M w /M N . Figure 3 show PDI vs t/τ 0 plot for all the three cases for different packing fractions. PDI is low at initial time due to presence of single NPs and very small clusters, it increases as time increases and reaches a maximum value at intermediate times. As time further increases, small aggregates start coalescing to form large aggregates which again leads to a decrease in PDI. At large times, PDI value saturates at PDI = 1 (corresponding to monodispersed case) for large η values which represent the presence of a single large percolated cluster. Also, at initial time PDI value is higher for larger η which indicates that aggregation is faster for high η. This further substantiates our point that initial time aggregation behavior is not captured for higher η values. Another important point to notice is that the time corresponding to the maxima of PDI (see Fig. 3) is directly proportional to the crossover time between two scaling regimes observed in aggregate size dynamics (see Fig. 1). Simulation starts with the initial state of PDI = 1, i.e. homogenous system, but as the system evolves PDI increases (see Fig. 3), which implies that the heterogeneity in the system increases. At the maxima of PDI, we have very heterogeneous system (many different size clusters with different charges) and we can predict the aggregation kinetics of this heterogeneous system by looking at the evolution of the system from this point onward. Now, comparing Fig. 1 and Fig. 3 we can easily conclude that time corresponding to the maxima of PDI is always greater than the crossover time between the first and second scaling regimes. Therefore, for a highly heterogeneous system, we can conclude that aggregation kinetics will have only one scaling regime.
Next, we plotted the radius of gyration square (R g 2 ) vs time in Fig. 4. Here we found an interesting non-monotonic behavior of R g 2 for high value of η and plateau like behavior in case of intermediate η. The peak arises for higher values of η due to formation of a single big aggregate that spans the whole simulation box, but we observe formation of fibril like structure at initial to intermediate times (see Fig. 2C). The system further evolves to minimize free energy, and therefore whole cluster evolves to compact form giving rise to a decrease in the R g 2 . For intermediate and small η values, system evolves very slowly as compared to high η values hence it gets sufficient time to restructure its clusters and thus we get a saturation or plateau like structure instead of maxima. Another interesting point to notice is that the peak in case of charge asymmetric is much broader as compared to charge symmetric cases. This tells us that the rearrangement of fibril like structure to a more compact structure will start very slowly in charge asymmetric case. This might be due to the strong Columbic attraction between the NPs in charge asymmetric case as compared to symmetric and size asymmetric case. Henceforth, this leads to high barriers in the potential energy landscape of charge asymmetric system and thus makes NPs rearrangement difficult and therefore slow.  www.nature.com/scientificreports www.nature.com/scientificreports/ As we can observe in Fig. 5A, aggregation is dominant and increasing the particle concentration increases the aggregation rate. Fluctuations in aggregate size, as clearly visible in Fig. 5A, might indicate the non-homogeneity of the mixture. Figure 5B shows the relative intensity of particles of different size at a particular time. We observed that peaks are shifting towards right that indicates the increase in the average size of the aggregates. Furthermore, intensity profile is getting broader also with time, which indicates that poly-dispersity of the system is increasing with time. At different times we also observed bimodal nature in intensity profile for example in Fig. 5B, see time = 30, 90 min. Bimodality confirms the non-homogeneity of the mixture. We did not observe two regimes in aggregation kinetics as observed in numerical simulation. One reason for this might be initial heterogeneity of the mixture, as explained earlier (i.e., starting from PDI maxima in simulations). Another possibility is that the NP used in the experiments are polymeric nanoparticles, which are coated with surfactant on the surface and they do not behave exactly like spherical NPs as we have assumed in numerical model. More controlled experiments are underway to study the system kinetics in more detail.

Conclusion
In conclusion, we studied the effect of charge and size asymmetry of NPs on their aggregation kinetics and compared it with the symmetric system. We found two aggregation regimes in the aggregation kinetics of three different systems viz. symmetric, charge asymmetric and charge asymmetric. The aggregation kinetics is strongly dependent on the concentration of the NPs in the system. At initial time for a particular concentration, Columbic forces between NPs decides the aggregation rate and therefore charge asymmetric system containing NPs of high valence aggregates faster than symmetric systems (both size symmetric and size asymmetric). This also leads to large decrease in charge/size ratio of newly formed aggregates in charge asymmetric system as compared to symmetric system and hence at intermediate times, aggregation become slower in asymmetric system. This happens because of the interplay between the rate of agglomeration of clusters and the size of agglomerating clusters. Further, NPs within a cluster rearrange themselves to minimize free energy and due to interplay between the time scale of aggregation and time scale of rearrangement leads to a non-monotonic behavior in R g 2 . In charge asymmetric case rearrangement of NPs is difficult within a cluster due to presence of strong Columbic interactions between oppositely charged NPs and therefore, R g 2 shows broader peak for charge asymmetric case as compared to symmetric case. We also performed few preliminary experiments on two sets of NPs and DLS studies confirm the dependence of particle concentration on aggregation kinetics. Due to organic group present on the nanoparticle surface, it is difficult to derive any power law from experiments and further experimental optimization is underway.