Dynamic self-assembly of photo-switchable nanoparticles

Abstract

Nanoparticles functionalized with photo-switchable ligands can be assembled into a broad range of structures by controlled light exposure. In particular, alternating light exposures provide the means to control formation of assemblies of various sizes and symmetries. Here, we use scaling arguments and Kinetic Monte Carlo simulations to study the evolution of reversible aggregates in a solution of periodically irradiated photo-switchable nanoparticles. Scaling estimates of the characteristic size and the mean separation of aggregates agree with the simulations. The transition probabilities in the Kinetic Monte Carlo scheme are derived from a renormalized master equation of the diffusion process. Simulations on a system of nanoparticles, interacting through Lennard-Jones pair potentials that change their character from repulsive to attractive depending on the light exposure, show that the slow diffusion of particles at low effective temperatures (where the attractions are much higher than the thermal energy) results in the formation of small, “kinetically frozen” aggregates. On the other hand, aggregation does not occur at high effective temperatures, where the attractions are comparable to the thermal energy. In the intermediate range of effective temperatures, “fluctuating” aggregates form that can be stabilized by applying light pulses of specific lengths and frequencies. The aggregate sizes increase by increasing the packing fraction and the aggregates undergo transition to a percolated “network” at high packing fractions. Light-control of inter-particle interactions can either inhibit or promote nucleation and growth, and can reduce gel and glass formation.

---

1 Introduction

Self-assembly of colloids and nanoparticles can generate a multitude of structures and patterns for a range of applications, including catalytic systems,^1,2 electronic devices^3,4 and colloidal inks.^5,6 One of the major complications in the field is the control of the assembly process such as to avoid kinetic trapping and defect formation.^7–9 Dynamic self-assembly^10–13 aims at producing ordered structures “on demand” by controlled energy influx into the system. Systems are driven out of their natural (equilibrium) states and forced to adapt ordered architectures, by changes in the chemical or field-based stimuli. This strategy, developed by Whitesides and Grzybowski^14,15 has been widely successful in the design of nanoscale components and systems, in cases where the conventional equilibrium design methods do not offer adequate and dynamic control over the formation of assemblies. Of the many possible candidates of the control variables, light stands out unique, as it is non-evasive and capable of instantaneous delivery at precise locations;^16,17 as opposed to the relatively slow chemical delivery methods for which diffusive limitations are often present. Recent experimental efforts on the light induced self-assembly (LISA) of photo-switchable materials¹⁸ have met with considerable success, and many applications such as in self-erasable papers,⁶catalytic switches,¹ and nanoscale temperature sensors¹⁹ are being considered.

One barrier to future developments in dynamic self-assembly is the lack of guiding general principles that would describe the behavior of systems at miniature length-scales and far-from-equilibrium. Despite long history of research—dating back to Onsager and Prigogine—non-equilibrium thermodynamics is a field still in its infancy.²⁰ In the absence of general description of far-from-equilibrium systems, one has to resort to computational approaches which, however, are system-specific and present a multitude of technical challenges associated with tracking the dynamics of infinitely many particles that are needed to explain the mesoscale behavior of the systems under study. A useful compromise is achieved by numerical schemes based on continuum models, coarse-grained particle models, and implicit solvent models. These methods enable the possibility of simulating larger systems for longer times, at the expense of some arbitrariness in the definition of time. Phase-field approaches,²¹ Coarse-grained Molecular Dynamics (MD),²² and Brownian dynamics methods²³ are examples of such schemes.

We propose here a scaling theory and a kinetic Monte Carlo (kMC) simulation scheme to study aggregation kinetics in a model non-equilibrium nanosystem of light-switchable nanoparticles (NPs) under periodic light forcing. As a case in point, we consider a variation of an experimental LISA process,¹⁶ where a solution of NPs functionalized with photo-switchable ligands is subjected to alternating exposures to UV and visible light (henceforth referred as “ON” and “OFF” states of light, respectively), leading to the conversion of associative NPs to non-associative NPs, and vice versa. Controlling the time duration of light exposure allows for a control of the sizes and the polydispersities of NP aggregates.⁶ At the same time, the periodic exposure is a test of the structural reversibility of these aggregates. This study of the phase behavior of “switchable” nanocolloids can thus provide a theoretical framework for studying “dynamic” materials whose properties derive from the ability of the components to change their interaction potentials in response to specific environmental stimuli. Also, the idea of periodic forcing of switchable NPs is motivated by the issues of energy efficiency in dynamic self-assembly where, as we show, achieving an ordered state does not necessarily require a constant flux of energy but only a properly timed sequences of “pulses” sufficient to maintain the ordered states.

2 Methodology

Fig. 1A illustrates a schematic of our model system, whereby we describe all relevant parameters of the model. The magnitude and the time span of light exposure are crucial in determining the resulting size and structure of the NP aggregates. In particular, in the case of no light exposure, all NPs are non-associative and the solution is akin to an ideal gas of particles. Under light irradiation and in the limit of very long aggregation times, all NPs are associative and one large aggregate is expected due to the Ostwald ripening effect. For intermediate values of the aggregation time, the system remains macroscopically homogeneous with local heterogeneities characterized by the separation between microscopic aggregates. Thus, it is reasonable to expect that the sizes of the NP aggregates and the separations between them can be controlled by varying the aggregation time τ_a and the disassembly time τ_d. In this spirit, the mean separation λ between the aggregates at the end of the aggregation phase, the mean free path of free NPs in the original solution λ₀, and the diameter of the NPs r₀, provide the relevant length scales for the problem (Fig. 1B). The packing fraction of NPs in the original solution η₀ is related to λ₀ by η₀ = (λ₀/r₀)⁻³. The choice of τ_a and τ_d must take into account the characteristic times for the non-associative (A) ⇔ associative (B) switching of individual NPs – these times are defined as τ_on and τ_off for the forward and backward transitions, respectively. In the scaling arguments that follow, we consider a scenario where the A ⇔ B transition is infinitely fast compared to assembly and disassembly times under different light exposures (i.e. τ_on ≪ τ_a and τ_off ≪ τ_d). We note that this assumption is realistic for the experimental LISA systems where the isomerization of the nanoparticles (e.g., ones covered with azobenzene-dithiol ligands¹⁶) under UV light is usually much faster than the assembly or disassembly of any sizable NP aggregate. Furthermore, τ_d is chosen such that the system is devoid of all heterogeneities and returns to a random distribution at the end of the disassembly phase. The latter criterion draws inspiration from the recent experiments on self-erasing materials,⁶ where complete structural reversibility was desired. As discussed later, we use the scaling estimates to guide the kMC simulations, where these simplifying assumptions are relaxed. In order to distinguish the scaling estimates from those used in simulations, we represent the scaling estimates of τ_d and τ_a by τ and τ′, respectively.

Fig. 1 Schematic of the LISA process with periodic light forcing. (A) Fraction of associative particles (green bold line) and light intensity (black dashed line) as a function of time. and τ_a are the time spent in the “disassembly” and the “aggregation” phase, where the light is turned “OFF” and “ON”, respectively. The associative ⇔ non-associative transition follows a first order kinetics, with the fraction of sticky particles given by the first order kinetics, f_B = 1 − e^−δt/τ_on (switch-on) and f_B = e^{−δt/τ_off} (switch-off), δt being the time spent since the last switch-on(off) event. (B) Characteristic length scales in the disassembly and the aggregation phase. In the disassembly phase, the length scale is given by the mean free path λ₀ = η^−1/3₀r₀ where η₀ is the packing fraction of particles in the original solution, and r₀ is the diameter of the particles. The aggregation phase has a characteristic length scale λ, which corresponds to the lattice parameter for ordered crystalline structures, and the mean distance between oscillating units for disordered fluctuating structures.

2.1 Scaling analysis

The characteristic time needed to completely destroy an NP aggregate is given by the familiar diffusion relation, τ = λ²/D, where D is the diffusion coefficient of the free NPs in solution. On the other hand, NP aggregation can be represented as a homogeneous nucleation process of the type B + B_N = B_N+1, where N(t) is the number of NPs in the aggregate at time t measured since the beginning of the aggregation phase. This description assumes active nucleation centers of number density n_c = 1/λ³. The number density of free NPs (B) can be given by the first order kinetics, with the rate constant given by the Smoluchowski relation, K(t) = 4πRD, where R ∼ r₀N^1/3 is the radius of the aggregate. The aggregation time can be estimated by using the life time of free particles obtained from the numerical solution of the kinetic equation dn(t)/dt = −K(t)n_cn(t) for the number density of free particles n(t) = n₀ − n_c, where n₀ = 1/λ³₀ is the number density of the particles in the original solution. The aggregation time thus obtained is given by†

Note that τ and τ′ are not independent of each other. For a given system, an increase in τ′ must be accompanied by an increase in τ by the same factor. Finally, the maximum number of particles in an aggregate, or “aggregation number”, can be estimated by the relation, N₀ = (λ/λ₀)³. As an illustration, we consider the typical range of experimental parameters:^1,16r₀ = 5–10 nm, D = 4 × 10⁻⁷ cm² s⁻¹, and n₀ = 10¹³–10¹⁴ cm⁻³. Then, one can estimate that the assembly and disassembly of aggregates comprising N₀ = 1 000 000 particles requires exposure times of τ′ = 250–5000 s and τ = 12–54 s, which is similar in magnitude to the experimentally determined values. Moreover, τ and τ′ scales as N^2/3₀ assuming that other system parameters are held fixed.

Naturally, the scaling estimates give only a rough approximation of the assembly/disassembly process and lack any information on the packing fraction and the structure of the aggregates, and the effects of the degree of particles' association and the diffusion rates. However, scaling estimates significantly reduce the microscopic parameter regime of a large number of model parameters, which needs to be explored to gain physical insights into the underlying process. Therefore, we use the scaling analysis to guide the kMC simulations, and to gain a detailed understanding of the assembly/disassembly process. Note that our current study is for a “forced dynamic self-assembly” process which does not possess a natural frequency and characteristic length scales of oscillations in the absence of periodic light switching, thus enabling us to develop analytical scaling estimates. In the more general case of “self-oscillating” systems with natural frequencies and characteristic length scales, preliminary simulations are required to establish the scaling laws.^24,25

2.2 Renormalization of diffusion equations

The traditional approach of solving diffusion problems uses either the numerical solution of the deterministic Smoluchowski diffusion equation or that of a related stochastic diffusion equation using Brownian dynamics. A common problem in both methods is that the choice of time-step is limited by the values of the interaction potential and temperature. This requirement is associated with the numerical stability of these methods, and if a proper timestep is not employed, one is often faced with the challenge of error propagation with simulation time. This problem is even more exacerbated if the interaction potentials are unbounded at a certain point in space due to some constraints (e.g., in the case of hard-core potentials), where the magnitude of time step must be small enough. Moreover, in a Brownian dynamics simulation, it is practically impossible to distinguish between the stochastic terms in the equations and those that appear due to numerical error. This accumulation of numerical error manifests in uncontrolled changes in the effective diffusion constant or the type of diffusion (normal to anomalous). Further, one particular realization of the stochastic “noise” term is often not relevant (and is physically unrealistic) and we are interested in the average quantities which are obtained by averaging over many trajectories. In the following, we discuss a renormalization approach that addresses all these problems. The key ingredient in this approach is that the equations of diffusion dynamics are renormalized to a master equation with transition rates identical to the Glauber definition of transition probabilities,²⁶ which can be efficiently solved by Monte Carlo simulations.

Let us consider the diffusion kinetics of a particle in one-dimension under the influence of an external force F(r). The probability density to find the particle at location r at time t, c(r,t), is given by the Smoluchowski equation,

∂_tc(r,t) = ∂_rj(r,t), (2)

where the flux j(r,t) is given as,

j(r,t) = D[∂_rc(r,t) − βF(r)c(r,t)], (3)

with β = 1/k_BT, where k_B and T are the Boltzmann constant and temperature, respectively. For numerical convenience, eqn (2) can be written as a difference equation of the form:where we assume space and time discretizations, r_m = am and t_n = nΔt. A central difference approximation is used for the space derivative of the flux and L_m is the linear operator for the time difference. In the Forward Euler scheme, L_mc=(cⁿ⁺¹_m − cⁿ_m)/Δt, where cⁿ_m indicate the probability density of the particle to be at location r_m at time t_n. Note that the simultaneous integral transformation and discretization as applied in eqn (4) ensures that the difference scheme is conservative meaning that the total probability stays constant.

Next, we introduce an auxiliary field,

w(r) = c(r)e^βU(r), (5)

where U(r) is the potential energy of particle at location r. Using F(r) = −dU(r)/dr in eqn (5), eqn (3) for the flux reduces to,

Integrating dw(r)/dr in the interval r ∈ [r_m,r_m+1], we get,

where we have used the fact that e^βU(r) in the above integrand changes much rapidly than the flux j(r) on the scale Δr = a. Using eqn (5) and eqn (7), we obtain the master equation,

Here, K(m + 1 → m) and K(m → m + 1) are the dimensionless transition rates defined as,

The resulting form of the flux term j_m+1/2 is a renormalized one, since the potential and temperature appear only inside the definition of coefficients, which are bounded. Thus, it is always possible to find a constant time step which is not a function of temperature or potential for which the numerical scheme will be stable (errors will not propagate in time). Further, it is obvious from eqn (12) and eqn (13) that the probability of transition into a non-physical zone, where U_m+1 → ∞, equals zero. Eqn (9) and eqn (10) for the dimensionless transition rates are identical with the definition of transition probability in Glauber dynamics.²⁶ Namely, for a change in state s_i → s_j with energies E_i and E_j, the transition probability is defined as,

Note that the transition rate defined in eqn (11) is symmetric and constrained in the interval [0,1]. At equilibrium (j_m+1/2 = 0) the flux equation (eqn (8)) reduces to the Boltzmann distribution, c_m+1 = c_mexp[-β(U_m+1 − U_m)].

2.3 Kinetic monte carlo (kMC) scheme

As derived above, the numerical solution of the diffusion equation through a difference scheme can be mapped to a finite step Monte Carlo scheme with step size equal to the spatial discretization of the difference scheme. However, to complete the transformation, we need a relation for the time step in terms of the number of Monte Carlo steps and the step size in space. As an example, we consider a one-dimensional diffusion problem with a constant drift force F, such that U(r) = −Fr. In a given Monte Carlo step, a particle may move from a position r_n to r_n+1 ∈ {r_n,r_n ± a} in a time duration [t_n,t_n+1 = t_n + Δt]. The transition probability for the movement is given according to eqn (11) as,where the last expression is valid for small F. Note that the probabilities defined above are for the cases when the direction of movement is first picked at random (with a probability 1/2 for both + and − directions), and the particle can stay at its original position with a probability 1/2. Thus, the mean drift is given by the relation:

Now, the mobility can be defined as,

Finally, using the Einstein's relation, D = μ/β, we can find the definition of timestep as, Δt = a²/4D. This result can be easily generalized to an arbitrary dimension d as,

The recipe of a kMC simulation using Glauber dynamics is similar to the conventional Metropolis scheme,²⁷ but with transition probabilities defined as

for a transition with the energy change ΔE. Further, the choice of particle displacements is limited to the continuum version of flip movements of lattice simulations. In other words, the attempted moves of particles are constrained to lie on a surface of a sphere of radius a centered at the “old” position, where a is the step size. We perform three-dimensional kMC simulations of a system of NPs, placed in a large cubic box with periodic boundary conditions. The simulation box is composed of cubic cells of size λ (Fig. 1B). The number of such cubic cells provides a rough measure of the number of aggregates expected in a simulation, as estimated by the scaling arguments. Note that the cubic grid is employed only to obtain better statistical measures and to minimize periodicity artifacts; it does not play any role in the evaluation of energetic interactions. We use a simulation box with 5³ = 125 cubic cells in all simulations. A kMC step consists of attempted moves of all particles, decided by the choice of the azimuthal angle ϑ = 2πu and the polar angle ϕ = acos(2v − 1),²⁸ where u and v are random variables in (1,3). The attempted moves are accepted with the transition probability given by eqn (16). The “new” positions of particles, if a move is accepted, is stored in a coordinate array, and is copied into the particle coordinates once all particles have moved. Each kMC step corresponds to a step in time given by eqn (15), that is, for d = 3,where τ₀ = r²₀/D is a convenient unit of time. We use a = 0.1r₀ in our simulations. This choice gives reasonable acceptance rates of particle movements. The degree of association is characterized by a dimensionless temperature θ ≡ T/T₀, where T and k_BT₀ are the temperature and the depth of the LJ potential, respectively. The interactions between particles at distance r are modeled using a shifted LJ potential of the form:

The tail-correction ε is chosen such that U(r_cut) = 0 for a cutoff distance r_cut. We choose r_cut = 2^1/6r₀ for interactions involving non-associative particles, to mimic purely repulsive hard cores of particles. A higher cutoff value of r_cut = 2.5r₀ is used for short-range attraction between associative particles, with the cutoff chosen for efficiency in the implementation.

The effective diffusion rates of NPs inside the aggregates are different from those of free NPs in dilute solution – characterized by the diffusion constant D. Indeed, the effective diffusion of NPs depends also on the local concentration, which is often described in a phenomenological manner by including arbitrary concentration dependence in D. Such ad hoc modification is not required in the current model, because the diffusion of NPs inside the aggregates is naturally reduced by the presence of nearby particles, as reflected in the acceptance rates of kMC moves. For example, at low temperatures, NPs diffuse slower inside the aggregates (low acceptance rates) than compared to the high temperature case or the free NPs in solution (high acceptance rates). Further, since our objective is to capture the true dynamics at short time-scales and not the equilibrium phase behavior at long time scales, we refrain from the use of collective moves that would necessitate the definition of multiple time scales.

2.4 Characterization of aggregates

The aggregate size distribution is sampled by counting the number of aggregates n_i, with an aggregation number i, which is defined as the number of particles in the aggregate (i = 1,2…). Two particles are considered to be the part of an aggregate if the distance between them is less than 1.5r₀.‡ Aggregate size distributions are used to compute the number average aggregate size and the weight average aggregate size . The polydispersity of the aggregates is defined as the ratio M_w/M_n and is equal to 1 for mono-dispersed aggregates. The radial distribution function is averaged over all particles in the simulation box at a given instant of time (not time-averaged) and obtained by the standard binning procedure.²⁹ Local bond order parameter q₆ for particle i is defined as,³⁰where the local bond orientational order parameter q̄_6m(i) is given by,

Here, Y_6m(r̂_ij) are the spherical harmonics of order six, r̂_ij are the unit vectors connecting the particle i with its nearest neighbors and N_b(i) is the total number of nearest neighbors of particle i. The nearest neighbors are found using a parameter-free neighbor search algorithm.³¹

3 Results and discussion

3.1 Aggregate size

To quantify the effects of the degree of association and nanoparticle concentration, we analyze the system behavior for a range of dimensionless temperatures, θ, and packing fractions of particles in the original solution, η₀ after one cycle of light exposure (Fig. 2 and Fig. 3). At low temperatures (θ = 0.05 in Fig. 2A), the system favors small, “kinetically frozen” aggregates with high degrees of symmetry. An increase in θ beyond this limit leads to an increase in size of the aggregates (θ = 0.2 in Fig. 2A) and the system favors “flexible” structures with fluctuating particles that remain in the aggregate. Increasing θ further results in “highly fluctuating” aggregates, which often dissociate and re-combine with the constituent NPs (θ = 0.6 in Fig. 2A). At even higher values of θ (θ = 1.0 in Fig. 2A), aggregation does not occur. This non-monotonic dependence on θ arises from the competition of energetic and kinetic effects. While the energetic interactions due to the association of particles decrease with an increase in θ, the effective diffusion rates increase with increasing θ. Therefore, the system tends to favor small, relatively mono-dispersed aggregates at low θ, which increase in size until θ reaches a threshold temperature, where the kinetic effects become dominant (between θ = 0.3 and θ = 0.4 in Fig. 3A). Beyond this threshold temperature, the aggregates begin to dissociate and re-combine and form highly fluctuating poly-dispersed structures, and no aggregation occurs at very high θ (θ > 0.7 in Fig. 3A). An increase in η₀ leads to an increase in the size of aggregates (Fig. 2B and Fig. 3B) eventually leading to the formation of percolated “network” structures (η₀ = 0.2963 in Fig. 2B).

Fig. 2 Typical simulation snapshots at the end of aggregation phase after one cycle of light exposure. Snapshots are zoomed in to show individual aggregates in detail and do not show the entire simulation box. The contrast (bold or faded) of particles is used to indicate the depth of particles from the top plane (close or farther). (A) For a range of θ with η₀ = 0.037. (B) For a range of η₀ with θ = 0.2. τ_d = τ, τ_a = τη^−1/3₀ and λ/λ₀ = 4 for all snapshots. The relaxation process is assumed to be infinitely fast.

Fig. 3 Number average aggregate size (M_n) and polydispersity (PD) of the aggregates at the end of aggregation phase after one cycle of light exposure. Variation of M_n and PD with (A) temperature θ ≡ T/T₀ for η₀ = 0.037, (B) with η₀ for θ = 0.2. The model parameters are τ_a = τη^−1/3₀, τ_d = τ, and λ/λ₀ = 4. The relaxation process is assumed to be infinitely fast. Error bars have roughly the size of the symbol and represent the numerical error as opposed to the error resulting from the presence of stochastic term in Brownian dynamics.

It is useful to note here that the effective temperature θ in this study is introduced as a measure of the degree of association and should not be interpreted as the actual temperature of the system at thermodynamic equilibrium. Indeed, the non-monotonic dependence of the aggregate size on θ is a consequence of the fact that the systems at different values of θ are subjected to equal exposure times, which are estimated by scaling arguments. We believe that the systems at low temperatures would tend to form much larger aggregates, if they are allowed sufficiently long time in the aggregation phase. However, the dissociation and re-combination of aggregates at high θ is expected to be present even at high temperatures in thermodynamic equilibrium, where the magnitude of inter-particle attractions is comparable to that of thermal fluctuations.

3.2 Aggregate structure

In Fig. 4A and Fig. 4B, we show the radial distribution function g(r/r₀) at the end of one cycle of light exposure. The sharp peaks in g(r/r₀) are indicative of the development of a local order inside the aggregates. The features in g(r/r₀) are enhanced with a moderate increase in θ (θ = 0.05–0.2 in Fig. 4A) and as η₀ increases (Fig. 4B) due to the increase in the number of layers of particles in the aggregates. However, the features disappear with a further increase in θ (θ = 0.4–0.8 in Fig. 4A) after the onset of dissociation and recombination of aggregates. It is difficult to quantify the crystallinity of aggregates of fewer than hundred particles observed in this study. The conventional procedure of crystal structure determination in simulations defines a particle to be “crystalline” or “solid-like” if the local neighborhood of the particle is identical (or strongly correlated) with its neighbors,^30,32 as determined by the correlations of the Steinhardt order arameters³³ or its modifications.³⁴ Such criterion is indeed not met in small nano-clusters, with few layers of particles and large surface to volume ratios. We use the distribution of local bond order parameter q₆ (Fig. 4C and Fig. 4D) as a probe to the development of crystalline order on an individual particle level.^30,32,33 We find that the distributions at low θ (Fig. 4C) and low η₀ (Fig. 4D) are shifted towards high q₆, indicating that the particles inside the aggregates are more crystalline in these cases. At low θ, particles pack tightly because of high inter-particle attractions, giving rise to a crystal-like ordering (or, for few particles, symmetric structures). On the other hand, at low η₀, the aggregates are less fluctuating as opposed to the high η₀ case where the crowding of particles leads to competing attractions from all the neighbors leading to higher flexibility. However, we did not observe a narrowing of q₆ distribution, as one would expect for perfect crystalline structures.³² With regards to the actual LISA experiments, we observe that a covalent cross-linking (of thiol groups) in addition to the short-range dipole–dipole attractions was necessary to induce perfect crystalline order,¹⁶ and high quality crystals were not formed in the absence of this crosslinking.⁶ However, the existence of such crosslinking renders the structures irreversible, as it is difficult to break the covalent bonds after it is formed.

Fig. 4 Radial Distribution function and local bond order parameter (q₆) distributions. Radial distribution function: (A) for a range of θ with η₀ = 0.037, (B) for a range of η₀ with θ = 0.2. Probability density distribution of q₆: (C) for a range of θ with η₀ = 0.037, (D) for a range of η₀ with θ = 0.2. The model parameters are τ_d = τ, τ_a = τη^−1/3₀ and λ/λ₀ = 4. The relaxation process is assumed to be infinitely fast.

3.3 Reversibility and relaxation

We have so far considered the final state of the system at the end of the first cycle of light exposure and assuming that the relaxation process is infinitely fast. Fig. 5 shows the system behavior after five cycles for different disassembly times, wherein the associative → non-associative transition occurs with a finite relaxation time such that the aggregates do not dissolve completely in the disassembly phase. The non-associative → associative transition is assumed to be infinitely fast, as in the previous discussion. This particular choice is guided by the experimental observation that the associative NPs in the aggregates take longer to switch to non-associative when the light is turned OFF, as opposed to the relatively fast switching of particles when the light is turned ON.^6,16 Aggregates evolve in a nonlinear way with larger rates of evolution at the beginning of the aggregation phase, and smaller rates towards the end of the aggregation phase. We observe that short disassembly times (τ_d = 0.05τ − 0.2τ) do not significantly affect the final aggregate size, when compared to the case where the light is kept ON for the same period (τ_d = 0). In fact, it may lead to an increase in the aggregate size, as observed for t > 1500τ₀ in the figure. This is because short disassembly times allow the aggregates to re-organize; they evolve at a faster rate when light is turned ON again. For intermediate disassembly times (τ_d = 0.5τ), the system evolves into smaller aggregates when compared to the τ_d = 0 case. For even higher disassembly times (τ_d = τ), the aggregation behavior is completely reversible and the aggregates evolve to approximately the same size after every cycle of light exposure. We performed this analysis for a range of θ values (not shown) and analyzed the temperature dependence of aggregate size at various values of τ_d for different cycles of light exposure. Typically, the choice of τ_d has smaller influence on aggregate size for high values of θ where diffusion effects are dominant meaning that the aggregates evolve to similar sizes after the light is switched off for a certain period, when compared to the case when the light is ON for the entire period. However, τ_d has a stronger influence for low values of θ, where the attractive forces are dominant, and an increase in τ_d beyond a certain value (τ_d > 0.2τ in Fig. 5) gives rise to smaller aggregates. Thus, the frequency of light exposure in combination with the exposure times can be used as a design variable to form aggregates of a desired size and reversibility.

Fig. 5 Aggregate size evolution with time and the effect of relaxation. The model parameters are η₀ = 0.037, τ_a = τη^−1/3₀, λ/λ₀ = 4, and θ = 0.2. The switch-on process is assumed to be infinitely fast but the switch-off process has a finite relaxation time τ_off = τ. The results are for five ON-OFF cycles.

The behavior shown in Fig. 5 exhibit striking similarities with the photo-initiated free-radical polymerization process of vinyl monomers.³⁵ The free-radical polymerization process contains an initiation step where the free radicals are formed from initiators (such as peroxides) followed by a propagation step wherein the polymer chains absorb these radicals to form longer chains. Polymerization terminates when these radicals combine leaving no free radicals. The rate of initiation and propagation steps is proportional to the intensity of light exposure and thus can be stopped by switching off light, while radicals may still combine thus reducing their concentration. The initiation/propagation and termination steps in this process are analogous to the aggregation and disassembly of NPs in the LISA process, except that in the free radical polymerization, the polymer chains once formed cannot degrade, as opposed to the forming and breaking of the assemblies in the LISA process. Therefore, the concentration of free radicals in solution shows a trend similar to that depicted in Fig. 5 for the aggregate size.

4 Conclusions

We develop here a combination of scaling theory and kinetic Monte Carlo scheme to analyze switchable self-assembly processes of nanoparticles. The model we implemented entailed simple potentials but its general features are representative of the experimental light-switchable nanosystems.^1,6,16 The system density and temperature (degree of association) are found to have a profound influence on the aggregation dynamics and on the resulting aggregate size distribution of the particles. Our model predicts the formation of fluctuating (or “flexible”) and “kinetically frozen” structures with varying degrees of association and packing fractions of the original solution, in agreement with the experiments.¹⁶ We also find network structures that resemble spinodal decomposition patterns for high packing fractions.³⁶ The competition between energetic and kinetic effects in light-controlled self-assembly affects local and global order of the forming structures. The local order of the aggregates, which is characterized here by the radial distribution function and the local bond parameter distribution, show a tendency to crystal-like ordering at low degrees of association and low packing fractions of the original solution. This provides a route to avoid frozen, disordered glass-like structures characteristic of many aggregating colloidal systems.⁷ Finally, we establish that the exposure times and the frequency of light exposure are important design variables to achieve structures with a desired size and reversibility.

The kinetic Monte Carlo scheme developed here provides several advantages over the conventional dynamical simulation methods (e.g., Brownian dynamics). Since particle displacements are accepted/rejected in a Monte Carlo fashion, as opposed to moving particles following a stochastic equation of motion, unphysical movements (e.g., violation of a hard core assumption) are not possible (these moves have zero acceptance). Further, the absence of a stochastic “noise” term resolves the computational difficulties associated with generating statistically independent trajectories with definitive mean properties. This advantage is particularly relevant for cases where the magnitude of the driving force is comparable to the stochastic noise in Brownian dynamics, such as in nucleation processes when Zeldovich factors are small or phase segregation process that occur close to critical points as well as the initial kinetics of phase segregation near spinodal decomposition lines. Finally, since the timestep is independent of the magnitude of the interaction forces, much longer time-steps can be employed than the conventional Brownian dynamics technique. It would be interesting to compare the accuracy and performance of this scheme for model systems to gain a detailed understanding of these differences.

The theoretical methods we described can be extended to study light induced aggregation in systems with more complex interactions, such as electrostatic interactions, providing a possibility to control re-entrant transitions characteristic of colloidal assembly^37,38 through the competition of short- and long- range interactions. In the future, it would be interesting to study systems of switchable colloids that can interact irreversibly (as in experiments where light-irradiation induces crosslinking of the NPs by covalent bonds). Our work shows that the kinetic Monte Carlo scheme for continuum (off-lattice) reaction-diffusion problems in soft-matter systems has great potential for large-scale simulations of non-equilibrium microstructures. The method can be extended to the design of specific morphologies and symmetries and to explore emergent phenomena by switching periodically and non-periodically the nature of the inter-particle interactions.

Acknowledgements

This work was supported by the Nonequilibrium Energy Research Center, which is an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0000989.

References

1. Y. H. Wei, S. B. Han, J. Kim, S. L. Soh and B. A. Grzybowski, J. Am. Chem. Soc., 2010, 132, 11018–11020.
2. T. H. Galow, U. Drechsler, J. A. Hanson and V. M. Rotello, Chem. Commun., 2002, 1076–1077.
3. A. J. Baca, J.-H. Ahn, Y. Sun, M. A. Meitl, E. Menard, H.-S. Kim, W. M. Choi, D.-H. Kim, Y. Huang and J. A. Rogers, Angew. Chem., Int. Ed., 2008, 47, 5524–5542.
4. D. V. Talapin, J. S. Lee, M. V. Kovalenko and E. V. Shevchenko, Chem. Rev., 2010, 110, 389–458.
5. J. C. Conrad, S. R. Ferreira, J. Yoshikawa, R. F. Shepherd, B. Y. Ahn and J. A. Lewis, Curr. Opin. Colloid Interface Sci., 2011, 16, 71–79.
6. R. Klajn, P. J. Wesson, K. J. M. Bishop and B. A. Grzybowski, Angew. Chem., Int. Ed., 2009, 48, 7035–7039.
7. M. E. Cates, Annales Henri Poincaré, 2003, 4, S647–S661.
8. K. A. Dawson, G. Foffi, F. Sciortino, P. Tartaglia and E. Zaccarelli, J. Phys.: Condens. Matter, 2001, 13, 9113–9126.
9. D. Chandler, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15111–15112.
10. B. A. Grzybowski and C. J. Campbell, Chem. Eng. Sci., 2004, 59, 1667–1676.
11. M. Fialkowski, K. J. M. Bishop, R. Klajn, S. K. Smoukov, C. J. Campbell and B. A. Grzybowski, J. Phys. Chem. B, 2006, 110, 2482–2496.
12. B. A. Grzybowski, C. E. Wilmer, J. Kim, K. P. Browne and K. J. M. Bishop, Soft Matter, 2009, 5, 1110–1128.
13. S. Mann, Nat. Mater., 2009, 8, 781–792.
14. G. M. Whitesides and B. Grzybowski, Science, 2002, 295, 2418–2421.
15. B. A. Grzybowski, H. A. Stone and G. M. Whitesides, Nature, 2000, 405, 1033–1036.
16. R. Klajn, K. J. M. Bishop and B. A. Grzybowski, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 10305–10309.
17. Q. Tran-Cong-Miyata, S. Nishigami, T. Ito, S. Komatsu and T. Norisuye, Nat. Mater., 2004, 3, 448–451.
18. M. M. Russew and S. Hecht, Adv. Mater., 2010, 22, 3348–3360.
19. R. Klajn, K. P. Browne, S. Soh and B. A. Grzybowski, Small, 2010, 6, 1385–1387.
20. M. Kröger, Phys. Rep.-Rev. Sec. Phys. Lett., 2004, 390, 453–551.
21. K. Thornton, J. Ågren and P. W. Voorhees, Acta Mater., 2003, 51, 5675–5710.
22. C. Reichhardt and C. J. O. Reichhardt, Phys. Rev. Lett., 2011, 106.
23. S. H. Northrup, S. A. Allison and J. A. Mccammon, J. Chem. Phys., 1984, 80, 1517–1526.
24. O. Kortlüke, V. N. Kuzovkov and W. von Niessen, Phys. Rev. Lett., 1999, 83, 3089–3092.
25. O. Kortlüke, V. N. Kuzovkov and W. von Niessen, Phys. Rev. E, 2002, 66.
26. R. J. Glauber, J. Math. Phys., 1963, 4, 294–307.
27. D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic Press, 2001.
28. E. W. Weisstein, From MathWorld–A Wolfram Web Resource.
29. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, USA, 1989.
30. P. R. ten Wolde, M. J. Ruiz-Montero and D. Frenkel, J. Chem. Phys., 1999, 110, 1591–1599.
31. J. A. van Meel, Dissertation, Van't Hoff Institute for Molecular Sciences (HIMS), 2009.
32. S. Auer and D. Frenkel, Adv. Polym. Sci., 2005, 173, 149–208.
33. P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Phys. Rev. B, 1983, 28, 784–805.
34. W. Lechner and C. Dellago, J. Chem. Phys., 2008, 129, 114707.
35. P. J. Flory, Principles of polymer chemistry, Cornell University Press, Ithaca, 1953.
36. P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino and D. A. Weitz, Nature, 2008, 453, 499–U494.
37. Y. L. Chen and K. S. Schweizer, J. Chem. Phys., 2004, 120, 7212–7222.
38. A. P. Hynninen and M. Dijkstra, Phys. Rev. Lett., 2005, 94.

---

Footnotes

† The numerical solution of the kinetic equation results in a prefactor of the order unity on the right hand side of eqn (1) that is not considered in the scaling estimates.

‡ The threshold distance must be chosen to be higher than the particle diameter r₀ and less than the cut-off distance of LJ interaction (r_cut = 2.5r₀ in our case). In general, the aggregates are loosely packed structures with mean distance between particles substantially higher than in the close packing limit (=r₀).

---