Generic model for tunable colloidal aggregation in multidirectional fields
Florian Kogler, Orlin D. Velev, Carol K. Hall, Sabine H. L. Klapp
GGeneric model for tunable colloidal aggregation in multidirectionalfields
Florian Kogler, ∗ a Orlin. D. Velev, b Carol K. Hall, b and Sabine H. L. Klapp a Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XXFirst published on the web Xth XXXXXXXXXX 200X
DOI: 10.1039/b000000x
Based on Brownian Dynamics computer simulations in two dimensions we investigate aggregation scenarios of colloidal particleswith directional interactions induced by multiple external fields. To this end we propose a model which allows continuous changein the particle interactions from point-dipole-like to patchy-like (with four patches). We show that, as a result of this change, thenon-equilibrium aggregation occurring at low densities and temperatures transforms from conventional diffusion-limited clusteraggregation (DLCA) to slippery DLCA involving rotating bonds; this is accompanied by a pronounced change of the underlyinglattice structure of the aggregates from square-like to hexagonal ordering. Increasing the temperature we find a transformation toa fluid phase, consistent with results of a simple mean-field density functional theory.
Recent progress in the synthesis and directional binding ofnanometer to micrometer sized patchy and anisotropic parti-cles makes possible the assembly of colloidal structures withmultiple directed bonds . The directional bonding can alsobe achieved by permanently embedded or field-induced dipoleand/or multipole moments allowing directional and selectiveparticle bonding . Within this class, particles with field-induced dipolar interactions are especially interesting be-cause switching the fields on and off is equivalent to switch-ing the particle interactions on and off. This means that ag-gregation mechanisms can be ’dialed in’. Furthermore,the orientation of inductive fields may be used to direct par-ticle aggregation . In consequence, such directed self-assembly processes may be exploited for the formation of newfunctional materials with specific and/or adjustable proper-ties. Hence, understanding the interplay between externallyinduced particle properties, external fields and thermodynamicconditions, e.g., temperature, is of fundamental interest inmodern material science, but also from a statistical physicspoint of view.An important subset of the many classes of self-assembledstructures are percolated colloidal networks, which are char-acterized by system-spanning cross-linked (patchy) parti-cle clusters that are realizable even at low volume frac-tions . Such network-like aggregates are consideredto be the underlying micro-structures of gels and have a Institute of Theoretical Physics, Technical University of Berlin, Harden-bergstr. 36, 10623 Berlin, Germany. E-mail: [email protected] b Department of Chemical and Biomolecular Engineering, North CarolinaState University, Raleigh, North Carolina, 27695-7905, USA. been intensively investigated in experiment and theory un-der equilibrium as well as non-equilibrium conditions .In the latter, qualitatively different aggregation mechanismscan be identified, namely diffusion limited cluster aggre-gation (DLCA) and reaction limited cluster aggregation(RLCA) . In the DLCA regime each particle collision leadsto the formation of a rigid and essentially (on the timescale ofthe experiment) unbreakable bond with fixed spatial orienta-tion. In contrast, in the RLCA regime the probability to forma rigid bond at collision is small. Systems with DLCA un-dergo irreversible dynamics and form fractal aggregates withspecific fractal dimensions D f ≈ . in continuous two-dimensional space . Such colloidal systems are consid-ered to be ’chemical gels’ and can be realized by having par-ticle interactions that are much stronger than k B T , preventingparticles from dissociating due to thermal fluctuations. Thisleads to a pronounced hindrance of structural reconfigurationof large particle aggregates . However, at higher tempera-tures these systems become ’physical gels’ where single par-ticles and larger substructures start to connect and disconnectfrequently. This strongly affects (increases) the fractal dimen-sion and finally allows the system to achieve its equilib-rium state.A recently introduced new type of DLCA, which accountsfor local rearrangements via flexible bonds, is slippery diffu-sion limited cluster aggregation (sDLCA) . Slippery bondsallow particles to move or rotate around each other as longas they stay in contact, meaning that bonds are still unbreak-able but can change their orientation. This additional degreeof freedom generates, at least in three-dimensional simula-tions , aggregates of the same fractal dimension as clas-sical DLCA but with a larger coordination number. ?? | a r X i v : . [ c ond - m a t . s o f t ] M a y LCA processes have been studied extensively in systemswith isotropically attractive particles but also in systemswith patchy particles bearing permanent and/or locally re-stricted interaction sites on their surfaces . In the latter,the spatial orientations of interaction sites can either be free torotate or fixed in space . When the orientationsof interaction sites are fixed in space, the associated ’chemicalgels’ undergo anisotropic diffusion limited aggregation whichyields a fractal dimension of D f ≈ . , lower than for theisotropic case. This situation occurs, e.g., due to the presenceof external fields or in lattice models , where motionis naturally restricted to certain directions.In the present paper we are particulary interested in theaggregation of colloids with field-induced multipolar inter-actions. Examples are capped (metal-coated dielectric parti-cles studied earlier by some of us ), where time-dependentelectric fields can induce quadrupolar-like interactions. Herewe consider even more complex interactions caused by crossed (orthogonal) fields. We briefly mention two examplesof possible experimental realizations of such systems. Thefirst one is a quasi two-dimensional system of suspended col-loidal particles, each composed of super-paramagnetic iron-oxide aggregates embedded in a polymer matrix, which hasbeen investigated experimentally by one of us . In this casecrossed external electric and magnetic fields, oriented in planebut perpendicular to each other, can be used to induce inde-pendent electric and magnetic dipole moments in the colloidsleading to a directed self-assembly process resulting in two-dimensional single-particle chain networks. A second possi-ble experimental and quasi two-dimensional system consistsof suspended colloidal particles under the influence of two in-plane orthogonal AC electric fields with a phase shift of π .The fields will polarize the particles’ ionic layer periodicallybut at different times due to their phase shift. By adjusting thefield frequencies and phases to the relevant timescales govern-ing particle diffusion and the relaxational dynamics of the po-larized ionic layer, two decoupled orthogonal dipole momentsin each particle can in principle be generated by this setup. Inboth cases, the crossed dipole moments might be characterizedas point-like or having a finite distance between their consti-tutive charges (or microscopic dipole moment distributions inthe magnetic case).Here, we investigate the structure formation in such systemsin a conceptional fashion by means of two-dimensional Brow-nian dynamics (BD) simulations of a generic particle model.The idea is to mimick externally-induced dipole moments viatwo pairs of screened Coulomb potentials that are decoupledto account either for magnetic and electric interactions or fortwo temporarily present electric interactions. The two chargesassociated with each pair are shifted outward from the par-ticle center, one parallel to the corresponding field and theother one anti-parallel. A sketch of such a particle with its Fig. 1 (Color online) Distribution of externally induced fictitious”charges” q inside a particle. Positions of charges are determined bythe vectors δ α k ∈ [ − δ e x , δ e x , − δ e y , δ e y ]) pointing either parallelor anti-parallel to the corresponding fields. internal arrangement of interaction centers is shown in Fig. 1.By changing the charge separation, we systematically investi-gate the (transient) structural ordering and aggregation behav-ior predicted by this model.Highlights of our results are the following: At very high in-teraction energies and large charge separations we find that theparticles undergo anisotropic diffusion limited cluster aggre-gation with rectangular local particle arrangements. Loweringthe charge separation shifts the model behavior to a slipperydiffusion limited aggregation (sDLCA) regime accompaniedby a sharp transition of the lattice structure from rectangu-lar to hexagonal. In the proximity of this transition we ob-serve long-lived or arrested frustrated structures consisting ofstrongly interconnected hexagonal and rectangular lattice do-mains connected with each other. We also show that, upon in-crease of the temperature, the systems enter a fluid state. Thecorresponding ’fluidization’ temperature turns out to be veryclose to the spinodal temperatures obtained from a mean-fielddensity functional theory.The rest of this paper is oranized as follows. In section 2we present our model and discuss relevant target quantities ob-tained in the BD simulations. Numerical results are describedin section 3, where we discuss first a specific low-temperature,low-density, state and then turn to the role of temperature anddensity. Finally, our conclusions are summarized in section 4. We consider a two-dimensional system of N soft spheres ofequal diameter σ . The soft sphere interactions are repulsive | ?? nd are modeled by a shifted and truncated (12,6) Lennard-Jones Potential U SS ( r ij ) = 4 (cid:15) (cid:0) ( σ/r ij ) − ( σ/r ij ) + 1 / (cid:1) (1)which is cut off at r c,SSij = 2 / σ . Here, r ij = | r j − r i | is theparticle center-to-center distance and (cid:15) sets the unit of energy.The crossed orthogonal external fields induce orthogonaldipole moments µ m = µ e m and µ e = µ e e which we termfor simplicity as ’magnetic’ and ’electric’ dipoles (althoughthe model is also appropriate for two electric moments). Thecoordinate frame is adjusted to coincide with the directionsof these moments so that e m = e x and e e = e y . In gen-eral these moments could have different absolute values butfor simplicity they are assumed to be equal. The two typesof dipole moments are also assumed to be independent fromeach other and interact only with dipole moments of the sametype on other particles. Intuitively, one would model the in-teraction energy between dipoles of particles 1 and 2 by thepoint-dipole potential U αdip ( r ) = µ α · µ α r − µ α · r )( µ α · r ) r , (2)where α indicates the dipole type as being either e or m . Dueto the constraint µ α (cid:107) µ α it follows that U αdip ( r ) = µ α µ α r (1 − r · e α ) r ) . (3)The resulting total dipolar interaction between two particles isthe sum of the dipolar potentials stemming from the magneticand electric dipoles, respectively. Using µ = | µ αi | and therelation ( r · e e + r · e m ) = r (which holds since µ e and µ m are orthogonal) we obtain U edip ( r ) + U mdip ( r ) = − µ r . (4)The resulting interaction on the right side of Eq. (4) is anisotropic, purely attractive interaction that lacks any kindof directional character. Therefore, the potential defined inEq. (4) can not generate any rectangular structures as observedin experiments . Moreover, having in mind that the field-induced spatial separation of charges in the ionic layer of asuspended particle is comparable to the particle diameter σ , apoint dipole model seems even more unreasonable.We therefore introduce an alternative model. Each dipolemoment µ α (with α = e, m ) is replaced by two oppositecharges − q α = q α which are shifted out of the particle cen-ter by a vector δ α k = ( − k δ e α , with k = 1 , . The vec-tor δ α k points either parallel ( k = 2 ) or antiparallel ( k = 1 )along the corresponding point dipole moment µ α . Indepen-dent of their shift or type, we set all charges to the same abso-lute value q = | q α k | = 2 . (cid:15)/σ ) − / . The choice of the value Fig. 2 (Color online) Normalized direction-dependent pair interac-tion U ( r ij ) [see Eq. (9)] between a particle in the center of the co-ordinate frame and a second particle (indicated as black circle in(a)) at various positions r ij for three different charge separations δ = 0 . , . , . σ corresponding to (a), (b), (c). Sterically ex-cluded areas are indicated by white circles. (d) Interaction energyat distance r ij = σ as function of φ , the angle measured in multiplesof π against the x-axis, for δ = 0 . σ (yellow), δ = 0 . σ (purple)and δ = 0 . σ (black). . is essentially arbitrary, as we will later normalize the in-teraction energy to eliminate the dependence of its magnitudeon the charge separation δ [see Eq. (8) below]. Charges k and l on different particles i and j interact via a Yukawa potential U α k α l ij ( r ij ) = − q exp ( − κr α k α l ij ) r α k α l ij (5)with r α k α l ij = | r j − r i + δ α l − δ α k | . The inverse screeninglength is choosen to κ = 4 . σ − and a radial cutoff r c = 4 . σ is applied. A schematic representation of the model with its in-ternal arrangement of ’charges’ is shown in Fig. 1. Due to theYukawa-like interaction our model lacks the long-range char-acter of true dipolar interactions. However, similar modelswith comparable interaction ranges have been used to describedipolar colloids in the context of discontinous molecular dy-namics simulations . We note that mimicking magneticdipoles via spatially separated ’charges’ appears somewhat ar-tificial from a physical point of view. However, here we couldthink of a particle with a strongly inhomogenous distributionof magnetic moments. Moreover, the idea behind our ansatz isnot to mimick a particular complex colloid, but rather to pro-vide a generic model for field-induced directional interactions.The arrangment of charges inside particles then results in a ?? | air-interaction U DIP ( r ij ) given by U DIP ( r ij ) = (cid:88) k,l =1 [ U e k e l ij ( r ij ) + U m k m l ij ( r ij )] . (6)In principle, U DIP ( r ij ) is a function of q and δ . To facilitatethe comparison between the interactions at different δ ( q ischosen to be constant), we normalize U DIP ( r ij ) according to ˜ U DIP ( r ij ) = U DIP ( r ij ) × u/U DIP ( σ e α ) (7)where the constant u = − . (cid:15) is calculated from the unnor-malized energy U DIP ( σ e α ) with model parameters δ = 0 . σ and q = 2 . (cid:15)/σ ) − / . This procedure ensures that the nor-malized energy between two particles at contact ( r ij = σ ) anddirection r ij = σ e α (pointing along one of the fields) has theconstant value u for all δ , that is ˜ U DIP ( σ e α ) = u. (8)The full pair interaction of our model is then given by U ( r ij ) = U SS ( r ij ) + ˜ U DIP ( r ij ) . (9)The resulting potential is illustrated in Fig. 2(a)-(c) for a parti-cle in the center of the coordinate frame and a second particleat various distances r ij and angles φ = arccos( r · e x /r ) with ’charge’ separations δ = 0 . , . , . σ . The value δ = 0 . σ is motivated by our simulation results presented inSec. 3.1. Sterically-excluded areas are shown in white and en-ergy values are color coded in units of (cid:15) . The weak anisotropyof the resulting particle interactions at small δ (where one es-sentially adds two dipolar potentials, see Eq. (4)) transformsto a patchy-like pattern by increasing δ . Energy minima be-come more and more locally restricted and particle interac-tions become directional in character. This is also seen inFig. 2(d) which gives the energy between two particles in con-tact as function of φ for different δ . From Fig. 2(d) we alsosee that, independent of the ’charge’ separation δ , the min-ima of the full interaction potential [see Eq. (9)] occur forconnection vectors r ij = σ e e and r ij = σ e m (i.e., pointingalong the fields). Note that this already holds for the unnor-malized energy given in. Eq. (6). Simulations are performedwith N = 1800 to particles at a range of reduced num-ber densities ρ ∗ = ρσ and temperatures T ∗ = k B T /(cid:15) , in asquare-shaped simulation cell with periodic boundary condi-tions. The equations of motion ˙ r i = T ∗ γ N (cid:88) j =1 ∇ U ( r ij ) + ζ i ( T ∗ ) (10)are solved via the Euler scheme with an integration stepwidth ∆ t = 10 − τ b , where τ b is the Brownian timescale defined by τ b = σ γ/T ∗ , γ = 1 . is a friction constant and ζ i ( T ∗ ) is aGaussian noise vector which acts on particle i and fulfills therelations (cid:104) ζ i (cid:105) = 0 and (cid:104) ζ i ( t ) ζ j ( t (cid:48) ) (cid:105) = 2 γk B T δ ij δ ( t − t (cid:48) ) .We perform simulations for up to τ b . To characterize the structure of the systems we consider sev-eral quantities. The first one is the mean coordination number ¯ z = 1 N N (cid:88) i =1 z i , (11)where z i is the number of neighbors of particle i and the sumis over all particles. In the following, two particles are consid-ered to be nearest neighbors if their center-to-center distanceis smaller than r b = 1 . σ .To identify local particle arrangements, the orientationalbond order parameter is of special importance. For particle k it is given by φ nk = 1 z k z k (cid:88) l =1 | exp( inθ klλ ) | (12)with z k being the number of neighbors and θ klλ = arccos( r kl · r kλ / ( r kl r kλ )) being the angle between the bond of particle k and its neighbor l measured against a randomly chosen bondof particle k to one of its neighbouring particles λ . Hence, φ nk = 0 for z k < . The integer value n determines the typeof order which is detected by this parameter. We concentrateon φ and φ to identify square (rectangular) and hexagonallattice types. Its ensemble average is calculated via Φ n = 1 N N (cid:88) i =1 φ ni . (13)The reversibility of ’bond’ formation and slipperyness of ex-isting bonds can be characterized by the bond and the bond-angle auto-correlation functions c b ( t ) and c a ( t ) . To evaluate c b ( t ) we assign a variable b ij ( t ) to each pair of particles ateach time step which is if the particles i and j are nearestneighbors or zero otherwise. The bond auto-correlation func-tion is then defined as c b ( t ) = (cid:104) b ij ( t ) b ij ( t ) (cid:105) , (14)where the brackets indicate an average over all pairs that arebonded at time t . The bond-angle auto-correlation function c a ( t ) is defined similarly by defining the unit vector a ij ( t ) = r ij ( t ) /r ij ( t ) , (15)such that c a ( t ) = (cid:104) − arccos ( a ij ( t ) · a ij ( t )) /π (cid:105) (16) | ?? ig. 3 (Color online) Simulation snapshots at ρ ∗ = 0 . and T ∗ = 0 . for (a) δ = 0 . σ , (b) δ = 0 . σ and (c) δ = 0 . σ . Particles arecolored according to their value of φ i . where we average again over all pairs. While c b gives theinformation on how stable bonds are over time, c a tells howstable their direction is over time. Note that in contrast to thetypical definition of correlation functions for stationary sys-tems , here the functions c b ( t ) and c a ( t ) are not independentof the time origin t .Finally, we consider the fractal dimension D f of parti-cle clusters, which is particularly important in the context ofDLCA. Clusters are defined as a set of particles with commonnext neighbors. The size of a cluster is then quantified by itsradius of gyration R g = 1 N cl N cl (cid:88) i =1 ( r i − ¯r ) , (17)where N cl is the number of particles in the cluster, and ¯r is thepostion of its center-of-mass. By plotting ln R g against ln N cl for different clusters, we extract the fractal dimension D f viathe relationship R g ∼ N /D f cl (see Ref. ). Our large-scale Brownian dynamics simulations show that thesystem is very sensitive to changes in temperature T ∗ , num-ber density ρ ∗ , and charge separation δ . In this large parameterspace we find a variety of different states ranging from smallfractal aggregates and single-chain structures at low tempera-tures to coarser, isolated or interconnected clusters at highertemperatures. In the following sections 3.1 - 3.3 we first dis-cuss the structure, the time correlation functions and the frac-tal dimensions at a low temperature and an intermediate den-sity, focussing on the impact of the model parameter δ . Insection 3.4 and 3.5 we then turn to the impact of temperatureand density. At first we study the system at low temperature T ∗ = 0 . andintermediate density ρ ∗ = 0 . for different charge separations δ . In Fig. 3 simulation snapshots for δ = 0 . σ, . σ, . σ at t = 300 τ b [see. Eq. (10) below] are shown, where τ b is theBrownian timescale. The colorcode reflects the orientationalbond order parameter φ i of each particle i . All three cases arecharacterized by clusters with irregular shapes. However, lo-cal particle arrangements differ strongly. While for δ = 0 . σ the particles aggregate in a hexagonal fashion, at δ = 0 . σ they aggregate into rectangular structures. At the intermedi-ate charge separation δ = 0 . σ , hexagonal order dominatesthe system; however, some clusters also reveal subsets of par-ticles in rectangular arrangements. A more quantitative de-scription is given by the orientational bond order parameters Φ shown in Fig. 4(a) as functions of δ . By increasing δ ,one observes a sharp transition at δ ≈ . σ from hexagonaltowards rectangular (square) order.The very presence of such a sharp transition can be ex-plained via energy arguments based on the δ -dependent pairpotential plotted in Figs. 2(a)-(d). To this end, we calculatethe energy U hexi ( δ ) = (cid:80) j =1 U ( r ij ) of a particle i with sixneighbors j , which are located in a hexagonal arrangement at’contact’ distance σ around i . Note that not all hexagonal con-figurations do have the same contact energy. This is due to theanisotropy of interactions, see Fig. 2(d). Therefore we con-sider a hexagonal configuration in which the contact energy isas low as possible (this configuration was found numerically).The dependence of this lowest contact energy U hexi ( δ ) on thecharge separation parameter is plotted in Fig. 5. Also shownis the corresponding energy U sqi ( δ ) = (cid:80) j =1 U ( r ij ) = 4 × u of a particle with four neighbors j located at distance σ in arectangular arrangement, i.e., in the energy minima around i ?? | ig. 4 (Color online) Results for simulations with N = 1800 attemperature T ∗ = 0 . and density ρ ∗ = 0 . . (a) Orientationalbond order parameters Φ for square (black) and Φ (yellow) forhexagonal particle arrangements. (b) Mean coordination number ¯ z as function of charge separation δ at times t = 100 , , τ b . (the quantity u was defined below Eq. (7)). Note that the en-ergy U sqi ( δ ) does not depend on δ according to Eq. (8). Asshown in Fig. 5, the two curves intersect at a ”critical” valueof δ = 0 . σ . Thus, the simple energy arguments alreadysuggest a transition between states with local hexagonal andsquare order, even though the predicted critical value is some-what larger than the value of δ = 0 . σ seen in the actualsimulations at finite temperature and density [see Fig. 4(a)].Further information is gained from the behavior of the meancoordination number as a function of δ plotted in Fig. 4(b) forthree different times t = 100 τ b , τ b and τ b . At all timesconsidered, ¯ z undergoes a steep decrease at δ ≈ . σ froma nearly constant value, ¯ z hex ≈ . , to a value ¯ z sq ≈ . .This behavior reflects, on the one hand, again the presenceof a sharp transition; on the other hand, the actual values of ¯ z hex (¯ z sq ) reveal the ”non-ideal” character of the aggregates interms of coordination numbers. For example, for δ > . σ we find that ¯ z and Φ decrease with δ , while Φ increases.However, this does not indicate a decline of the rectangu-lar order; it rather results from an increasing amount of par-ticles residing in chains oriented either in x- or y-direction.The coordination number z i of a particle i in such a chain is ≤ , leading to a mean coordination number ¯ z < . Further-more, the parameters φ i and φ i [see Eq. (12)] become unity Fig. 5 (Color online) Minimum energy of a particle with six neigh-bors in hexagonal arrangement as function of δ (black) and energyfor a particle in rectangular arrangement with 4 neighbors (red). for a particle forming exactly two bonds under an angle of π (straight chain). This does not affect Φ , which is alreadylarge at δ > . σ , but significantly increases Φ . Finally,the counter-intuitive decrease of Φ with δ results from theincreasing amount of particles with only one neighbor (e.g.,ends of chains appearing white in Fig. 3(c)). These particlesyield no contribution to Φ [see Eq. (12)].The ”non-ideal” values of ¯ z hex and ¯ z sq also explain whyour energy argument for the location of the hexagonal-to-square transition, which was based on ideal arrangements withsix and four neighbors, respectively, does yield the transitionvalue δ = 0 . σ rather than δ = 0 . σ obtained from sim-ulation. We can now reformulate the argument by using theactual mean coordination numbers extracted from our simu-lations, ¯ z sq = 3 . (instead of ) and ¯ z hex = 4 . (insteadof ). Following the calculations for the ideal arrangementsdescribed before, the energy of the square-like arrangement is U sq = 3 . × u . For the hexagonal arrangement, we use the av-erage minimum energy with either z i = 4 or z i = 5 neighbors,yielding ¯ U h (¯ z hex , δ ) = ( U hex (4 , δ ) + U hex (5 , δ )) / . The re-sulting critical value of the charge separation is δ ≈ . σ ,which coincides nicely with the transition value observed inour simulations. Although the local structures characterized by ¯ z and Φ per-sist, in general, over the simulation times considered, we arestill facing a transient (out-of-equilibrium) structure forma-tion as seen, e.g., from the slight increase of ¯ z with time inFig. 4(b). This raises a question about the typical ”lifetime” ofthe aggregates.To this end we now consider dynamical properties, namelythe bond and bond-angle auto-correlation functions, c b ( t ) and c a ( t ) . It is not reasonable to extract decay rates from thesefunctions (as it is usually done) because in transient states, de-cay rates are, strictly speaking, functions of time themselves.Still, it is interesting to see whether the temporal correlation | ?? ig. 6 (Color online) Time correlation functions obtained from sim-ulations with N = 1800 at temperature T ∗ = 0 . and density ρ ∗ = 0 . . (a) [(b)] Time evolution of the bond [angle] autocor-relation function c b ( t ) [ c a ( t ) ] for three different charge separations δ = 0 . σ, . σ and . σ colored in yellow, purple and black re-spectively. of bonds (bond angles) for different δ allows us to distinguishbetween qualitatively different aggregation regimes.Numerical results for c b ( t ) and c a ( t ) are plotted inFigs. 6(a) and (b), respectively, where we consider a large timerange up to t ≈ τ b . The time axis starts at the finite timewhen all the systems have formed stable aggregates. The datain Figs. 6(a) and (b) pertain to three representative values ofthe charge separation parameter related to the hexagonal struc-tures ( δ = 0 . σ ), rectangular structures ( δ = 0 . σ ), and to thetransition region ( δ = 0 . σ ). In the square regime ( δ = 0 . σ )the decay of both c b ( t ) and c a ( t ) is almost identical and veryslow. From this we conclude that the square regime is charac-terized by almost unbreakable bonds with fixed orientations.This is different in the hexagonal regime ( δ = 0 . σ ) where c b ( t ) remains nearly constant even after long times (meaningthat bond-breaking is very unlikely), while c a ( t ) decays muchfaster. Thus, the direction of bonds are less restricted. Weinterprete this behavior as evidence that two particles, thoughbeing bonded, are still able to rotate around each other to someextent. This is a characteristic feature of slippery bonds. Fi-nally, in the transition regime ( δ = 0 . σ ) both functions c b ( t ) and c a ( t ) decay significantly faster than in the other cases,with the decay of the bond-angle correlation function beingeven more pronounced. In that sense we may consider thebonds in the transition region also as slippery (although lesslong-lived than in the other cases).We conclude that the different structural regimes identifiedin the preceding section are indeed characterized by differentrelaxational dynamics. Moreover, all of the observed aggre-gates have lifetimes of at least several hundered τ b . Such long-lived bonds are indicative of diffusion limited cluster-clusteraggregation. In the next section we therefore consider the frac-tal dimension. Fig. 7 (Color online) Fractal dimension D f as a function of chargeseparation at ρ ∗ = 0 . and T ∗ = 0 . . At δ = 0 . σ we find abimodal distribution of fractal dimension with peaks at D f = 1 . (solid line) and . (dashed line). In Fig. 7 the fractal dimension D f is shown as a function of δ at time t = 250 τ b , density ρ ∗ = 0 . and temperature T ∗ =0 . . We find that D f increases slightly with δ but remainsin a range between . and . , except at δ = 0 . σ . There,the fractal dimension exhibits a bimodal distribution, takingvalues between D f ≈ . and D f ≈ . (dashed line inFig. 7).Despite these variations, the values of D f found here aresignificantly smaller than the fractal dimension D f = 1 . ob-served in earlier studies of DLCA in two-dimensional contin-uous (off-lattice) systems . Except for the case δ = 0 . σ ,the values in Fig. 7 are comparable with previous findings forDLCA in two-dimensional lattice systems and systems withspatial or interaction anisotropies . The present system isindeed anisotropic in the sense that the external fields imposepreferences on the directions of particle bonds and thereforealso on the orientations of aggregates. This effect is mostpronounced in the rectangular regime ( δ = 0 . σ ). There-fore, it is plausible that our system undergoes a special caseof anisotropic DLCA, in (quantitative) accordance with ex-perimental results and theoretical predictions . Weshould note that, due to our simulation method, the clustersizes (typically involving − particles) are relativelysmall compared to the particle numbers considered in the lit-erature ( particles) . A more detailed analysis of theimpact of anisotropic interactions on the fractal structure isbeyond the scope of this study.We also relate our findings to the newer concept of slip-pery DLCA , where the bonds are essentially unbreakablebut able to rotate. Indeed, as discussed in section 3.2, bondsare slippery in nature for small δ in the hexagonal regime.For three-dimensional systems it has been reported thatthe fractal dimension D f remains the same for slippery andclassical DLCA, while the mean coordination number ¯ z dif-fers. Specifically, ¯ z is significantly higher for sDLCA . ?? | ig. 8 (Color online) Temperature dependence of the system properties at density ρ ∗ = 0 . for charge separations δ = 0 . , . , . σ coloredin yellow, purple and black, respectively. (a) Fractal dimension D f evaluated at t ≈ τ b , (b) Mean coordination number, (c) Orientationalorder parameters Φ and Φ . The same observation emerges when we consider our valuesof ¯ z plotted in Fig. 4(b), from which one sees a pronounceddecrease of ¯ z upon entering the square (DLCA) regime. How-ever, in contrast to earlier studies we find D f to slightly in-crease with δ , especially in the hexagonal regime. We in-terpret this behavior as a consequence of the fact that bind-ing energies in the hexagonal regime decrease with increasingvalues of δ , while they remain constant in the square regime(see Fig. 5). The corresponding stability of bonds should becorrelated to the binding energies which explains the slightlyincreasing values of D f in the hexagonal regime.Finally, in the transition region ( δ = 0 . σ ) we found a bi-modal distribution of the fractal dimensions D f with maximaat D f ≈ . and D f ≈ . . This second maximum corre-sponds to only ≈ of the considered cases (twelve inde-pendent simulation runs). The first maximum at D f ≈ . therefore clearly dominates and fits nicely to the functionaldependence of D f on δ (see Fig. 7). We assume that theless frequent peak results from a switching of the local struc-tures between hexagonal and rectangular arrangements, whichis accompanied by a significantly larger bond-breaking prob-ability (see Fig. 6(c)). Again this allows compactification ofaggregates and increases the fractal dimension in the transitionregime. Diffusion limited aggregation is restricted to systems with at-tractive particle interactions much stronger than k B T . By in-creasing the temperature sufficiently, thermal fluctuations be-come able to break bonds which results in a faster decay of thethe bond auto-correlation functions and a compactification ofaggregates. Indeed, for square lattice models it was found that D f is a monotonically increasing function of temperature .In Fig. 8(a) the fractal dimension D f of the present model isplotted as a function of temperature T ∗ for charge separations δ = 0 . σ, . σ and . σ . We first concentrate on the case δ = 0 . σ , correspondingto the square regime at low T ∗ . In the range of very low tem-peratures T ∗ < . , the fractal dimension is small and staysessentially constant. Increasing T ∗ towards slightly larger val-ues then leads to an increase of D f , reflecting the (expected)compactification. This increase of D f is accompanied by anincrease of the mean coordination number ¯ z [see Fig. 8(b)]within the temperature range considered, indicating the grow-ing number of bonds due to local and global structural re-configurations. The corresponding changes in the stability ofthe bonds are illustrated in Fig. 9, where we have plotted thetime evolution of c b ( t ) for several temperatures (at δ = 0 . σ ).Clearly, the decay of c b ( t ) becomes faster for higher tempera-tures. This is the reason why structural reconfigurations and,in consequence, compactification of aggregates becomes pos-sible.These trends persist until T ∗ f,sq ≈ . , beyond which thesystem at δ = 0 . σ starts to behave in a qualitatively differ-ent way. The mean coodination number ¯ z displays a maxi-mum and subsequently a rapid decay. We also find that thefractal dimension has not yet reached its maximum value at T ∗ f,sq ; this maximum occurs at the slightly larger temperature T ∗ ≈ . (see Fig. 8(a)). This ’delay’ of D f can be under- Fig. 9 (Color online) Bond auto correlation function c b ( t ) for dif-ferent temperatures T ∗ at charge separation δ = 0 . σ and density ρ ∗ = 0 . . | ?? ig. 10 (Color online) Simulation snapshots at ρ ∗ = 0 . with δ = 0 . σ at (a) T ∗ = 0 . , (b) T ∗ = 0 . , and (c) T ∗ = 0 . . Particles arecolored according to their value φ i . stood from the fact that, upon the entrance of bond-breaking,filigree parts of the aggregates are more likely affected thanmore compact ones. Hence, the fraction of ’compact’ smallaggregates still grows. Even more important, the function Φ ( T ∗ ) in Fig. 8(c) displays a pronounced decay of rectan-gular order for T ∗ > T f,sq . From the sum of these indicationswe conclude that, at temperatures higher than T ∗ f,sq ≈ . ,the system transforms into a (stable or metastable) fluid phase.In this fluid phase, the overall structure starts to become ho-mogeneous and isotropic, while the local structures involveonly a small number of bonds with short bond-life times.For the system at δ = 0 . σ (hexagonal structure at low T ∗ ),an estimate of the ”fluidization” temperature T ∗ f,hex based onthe behavior of order parameters, coordination number andfractal dimension is more speculative. Nevertheless, the datasuggest that T ∗ f,hex > T ∗ f,sq . This is indicated, first, by thefact that Φ ( T ∗ ) decays only very slowly with temperatureuntil T ∗ ≈ . (see Fig. 8(c)). Second, the mean coordinationnumber shows only a weak maximum (and no fast decay after-wards) compared to the case δ = 0 . σ . Third, the fractal di-mension keeps increasing with T ∗ for all considered tempera-tures T ∗ < . . Therefore we conclude that T ∗ f,hex > . . Weunderstand this higher fluidization temperature at δ = 0 . σ from the fact that binding energies in hexagonal structures arelarger; therefore, higher coupling energies must be overcome.To further justify these interpretations, particularly theemergence of fluid phases, we performed a stability analy-sis of the homogenous isotropic high temperature state basedon mean-field density functional theory (DFT). Specifically,we consider the isothermal compressibility χ T . Positivevalues of χ T imply that the homogeneous (fluid) phase isstable, whereas negative values indicate that this phase isunstable. Specifically, the instability arises against long-wavelength density fluctuations, i.e. condensation. According to Kirkwood-Buff theory one has χ − T ∝ − ρ ˜ c ( k = 0) , (18)where ˜ c (0) is the Fourier transform of the direct correlationfunction (DCF) c ( r ) in the limit of long-wavelengths ( k → ). We approximate the DCF for distances r ij > σ accordingto a mean field (MF) approximation, that is c MF ( r ) = − ( k B T ) − U ( r ) , r > σ, (19)and use the Percus-Yevick DCF c HS ( r ) of a pure hard-sphere fluid for | r | ≤ σ . The full DCF is then given by c ( r ) = c HS ( r ) + c MF ( r ) . (20)In Fig. 11 we present numerical results for the expression − ρ ˜ c (0) at ρ ∗ = 0 . as function of temperature. At low T ∗ ,all systems are characterized by negative values of − ρ ˜ c (0) .This indicates that the homogeneous isotropic phase is unsta-ble, consistent with the results of our simulations. Upon in-creasing T ∗ the mean-field compressibility χ T then becomes Fig. 11 (Color online) Numerical solutions to Eq.18 as function of T ∗ for density ρ ∗ = 0 . and charge separations δ = 0 . , . , . σ colored in yellow, purple and black, respectively. ?? | ndeed positive for all charge separations considered. Specif-ically, for δ = 0 . σ the change of sign (related to a ”spinodalpoint”) occurs at T ∗ f,sq = 0 . and for δ = 0 . σ at the muchhigher temperature T ∗ f,hex = 0 . . These values are in surpris-ingly good agreement with our estimates for the ”fluidization”temperatures based on the order parameter plots.The case δ = 0 . σ is again different. Here we find[see Fig. 8(b)] that, starting from low temperatures inside theDLCA regime, the mean coordination number monotonicallydecreases. However, this does not indicate ”fluidization” butrather a gradual transition from a state with dominant hexago-nal order towards a mixed state comprised of coexisting clus-ters with local hexagonal and square-like order. Indeed, [seeFig. 8(c)], the orientational order parameters Φ and Φ revealthat the fraction of particles bound in square clusters increaseswith T ∗ and finally overtakes the fraction of particles involvedin hexagonal clusters at T ∗ ≈ . . Corresponding snapshotsof simulation results are shown in Fig. 10. At all temperaturesconsidered one observes separated clusters. With increasingtemperature their shape becomes more regular, while the lo-cal rectangular order becomes more pronounced. Finally, at T ∗ = 0 . the fractal dimension D f and the square orderparameter Φ reach their maximum values, suggesting a ”flu-idization” similar to the behavior observed at other values of δ .Interestingly, our stability analysis [see Eq. (18)] indicates aninstability at the same temperature T ∗ f = 0 . . With this sur-prisingly accurate agreement between theory and simulation,we conclude that in the transition regime ( δ = 0 . σ ), increas-ing thermal fluctuations first push the system from a domi-nantly hexagonal state into a rectangular one, which then en-ters a metastable fluid phase after passing the ”spinodal point”. In this section we revisit the system behavior at the low tem-perature T ∗ = 0 . , but consider different densities in therange ρ ∗ ≤ . . Whereas low-density systems at T ∗ = 0 . display DLCA as discussed in sections 3 A-C, this aggrega-tion mechanism is expected to disappear at higher densities:here, the particles are just unable to diffuse sufficiently freely.Rather, the particles will very frequently collide and then im-mediately form rigid bonds. A typical structure at the high-est density considered, ρ ∗ = 0 . , and separation parameter δ = 0 . σ is shown in Fig. 12. Clearly, the system is per-colated, that is, the particles form a single, system-spanningcluster. Interestingly, this cluster is composed of extended re-gions characterized by either square-like order or hexagonalorder. We note that, at δ = 0 . σ , simultaneous appearanceof clusters with both types of order also occurs at low den-sities and higher temperatures (see section 3.4). However, atthe high density considered here the regions of each type arelarger and the particle arrangements are much more regular Fig. 12 (Color online) System at T ∗ = 0 . , density ρ ∗ = 0 . and δ = 0 . σ . The colorcode gives the orientational bond-orderparameter φ i of each particle i . (i.e., there are less defects).To better understand the impact of the density on the clusterstructures we plot in Fig. 13(a) the orientational bond orderparameters Φ and Φ as functions of ρ ∗ for δ = 0 . σ (at δ = 0 . σ and δ = 0 . σ the order parameters are essentiallyindependent of the density). From Fig. 13(a) it is seen thatthe amount of rectangular (hexagonal) order sharply increases(decreases) at a density of ρ ∗ ≈ . . This is a surprisingresult as one would expect that, upon compressing the system,close-packed, hexagonal structures rather become more likely.However, at the low temperature considered here, structuralreorganization is strongly hindered.We also note that all of the systems investigated at densi-ties ρ ∗ > . turned out to be percolated (suggesting that thevalue ρ ∗ = 0 . is indeed related to the percolation transi-tion). It thus seems that the percolation tends to stabilize theinitially formed square-lattice symmetry, as the subsequent re-organization is hindered by the lack of mobility. In effect, weare faced with quenched states that could not densify withinthe time domain studied. This interpretation is also consistentwith the decrease of the mean coordination number once thesystem is percolated ( ρ ∗ > . ) as shown in Fig. 13(b). In this work we propose a new model for field-directed ag-gregation of colloidal particles with anisotropic interactionsinduced by external fields. The model was inspired by recentexperimental work on novel colloidal particles inwhich external fields can induce two, essentially decoupled,dipoles. |1–
Recent progress in the synthesis and directional binding ofnanometer to micrometer sized patchy and anisotropic parti-cles makes possible the assembly of colloidal structures withmultiple directed bonds . The directional bonding can alsobe achieved by permanently embedded or field-induced dipoleand/or multipole moments allowing directional and selectiveparticle bonding . Within this class, particles with field-induced dipolar interactions are especially interesting be-cause switching the fields on and off is equivalent to switch-ing the particle interactions on and off. This means that ag-gregation mechanisms can be ’dialed in’. Furthermore,the orientation of inductive fields may be used to direct par-ticle aggregation . In consequence, such directed self-assembly processes may be exploited for the formation of newfunctional materials with specific and/or adjustable proper-ties. Hence, understanding the interplay between externallyinduced particle properties, external fields and thermodynamicconditions, e.g., temperature, is of fundamental interest inmodern material science, but also from a statistical physicspoint of view.An important subset of the many classes of self-assembledstructures are percolated colloidal networks, which are char-acterized by system-spanning cross-linked (patchy) parti-cle clusters that are realizable even at low volume frac-tions . Such network-like aggregates are consideredto be the underlying micro-structures of gels and have a Institute of Theoretical Physics, Technical University of Berlin, Harden-bergstr. 36, 10623 Berlin, Germany. E-mail: [email protected] b Department of Chemical and Biomolecular Engineering, North CarolinaState University, Raleigh, North Carolina, 27695-7905, USA. been intensively investigated in experiment and theory un-der equilibrium as well as non-equilibrium conditions .In the latter, qualitatively different aggregation mechanismscan be identified, namely diffusion limited cluster aggre-gation (DLCA) and reaction limited cluster aggregation(RLCA) . In the DLCA regime each particle collision leadsto the formation of a rigid and essentially (on the timescale ofthe experiment) unbreakable bond with fixed spatial orienta-tion. In contrast, in the RLCA regime the probability to forma rigid bond at collision is small. Systems with DLCA un-dergo irreversible dynamics and form fractal aggregates withspecific fractal dimensions D f ≈ . in continuous two-dimensional space . Such colloidal systems are consid-ered to be ’chemical gels’ and can be realized by having par-ticle interactions that are much stronger than k B T , preventingparticles from dissociating due to thermal fluctuations. Thisleads to a pronounced hindrance of structural reconfigurationof large particle aggregates . However, at higher tempera-tures these systems become ’physical gels’ where single par-ticles and larger substructures start to connect and disconnectfrequently. This strongly affects (increases) the fractal dimen-sion and finally allows the system to achieve its equilib-rium state.A recently introduced new type of DLCA, which accountsfor local rearrangements via flexible bonds, is slippery diffu-sion limited cluster aggregation (sDLCA) . Slippery bondsallow particles to move or rotate around each other as longas they stay in contact, meaning that bonds are still unbreak-able but can change their orientation. This additional degreeof freedom generates, at least in three-dimensional simula-tions , aggregates of the same fractal dimension as clas-sical DLCA but with a larger coordination number. ?? | a r X i v : . [ c ond - m a t . s o f t ] M a y LCA processes have been studied extensively in systemswith isotropically attractive particles but also in systemswith patchy particles bearing permanent and/or locally re-stricted interaction sites on their surfaces . In the latter,the spatial orientations of interaction sites can either be free torotate or fixed in space . When the orientationsof interaction sites are fixed in space, the associated ’chemicalgels’ undergo anisotropic diffusion limited aggregation whichyields a fractal dimension of D f ≈ . , lower than for theisotropic case. This situation occurs, e.g., due to the presenceof external fields or in lattice models , where motionis naturally restricted to certain directions.In the present paper we are particulary interested in theaggregation of colloids with field-induced multipolar inter-actions. Examples are capped (metal-coated dielectric parti-cles studied earlier by some of us ), where time-dependentelectric fields can induce quadrupolar-like interactions. Herewe consider even more complex interactions caused by crossed (orthogonal) fields. We briefly mention two examplesof possible experimental realizations of such systems. Thefirst one is a quasi two-dimensional system of suspended col-loidal particles, each composed of super-paramagnetic iron-oxide aggregates embedded in a polymer matrix, which hasbeen investigated experimentally by one of us . In this casecrossed external electric and magnetic fields, oriented in planebut perpendicular to each other, can be used to induce inde-pendent electric and magnetic dipole moments in the colloidsleading to a directed self-assembly process resulting in two-dimensional single-particle chain networks. A second possi-ble experimental and quasi two-dimensional system consistsof suspended colloidal particles under the influence of two in-plane orthogonal AC electric fields with a phase shift of π .The fields will polarize the particles’ ionic layer periodicallybut at different times due to their phase shift. By adjusting thefield frequencies and phases to the relevant timescales govern-ing particle diffusion and the relaxational dynamics of the po-larized ionic layer, two decoupled orthogonal dipole momentsin each particle can in principle be generated by this setup. Inboth cases, the crossed dipole moments might be characterizedas point-like or having a finite distance between their consti-tutive charges (or microscopic dipole moment distributions inthe magnetic case).Here, we investigate the structure formation in such systemsin a conceptional fashion by means of two-dimensional Brow-nian dynamics (BD) simulations of a generic particle model.The idea is to mimick externally-induced dipole moments viatwo pairs of screened Coulomb potentials that are decoupledto account either for magnetic and electric interactions or fortwo temporarily present electric interactions. The two chargesassociated with each pair are shifted outward from the par-ticle center, one parallel to the corresponding field and theother one anti-parallel. A sketch of such a particle with its Fig. 1 (Color online) Distribution of externally induced fictitious”charges” q inside a particle. Positions of charges are determined bythe vectors δ α k ∈ [ − δ e x , δ e x , − δ e y , δ e y ]) pointing either parallelor anti-parallel to the corresponding fields. internal arrangement of interaction centers is shown in Fig. 1.By changing the charge separation, we systematically investi-gate the (transient) structural ordering and aggregation behav-ior predicted by this model.Highlights of our results are the following: At very high in-teraction energies and large charge separations we find that theparticles undergo anisotropic diffusion limited cluster aggre-gation with rectangular local particle arrangements. Loweringthe charge separation shifts the model behavior to a slipperydiffusion limited aggregation (sDLCA) regime accompaniedby a sharp transition of the lattice structure from rectangu-lar to hexagonal. In the proximity of this transition we ob-serve long-lived or arrested frustrated structures consisting ofstrongly interconnected hexagonal and rectangular lattice do-mains connected with each other. We also show that, upon in-crease of the temperature, the systems enter a fluid state. Thecorresponding ’fluidization’ temperature turns out to be veryclose to the spinodal temperatures obtained from a mean-fielddensity functional theory.The rest of this paper is oranized as follows. In section 2we present our model and discuss relevant target quantities ob-tained in the BD simulations. Numerical results are describedin section 3, where we discuss first a specific low-temperature,low-density, state and then turn to the role of temperature anddensity. Finally, our conclusions are summarized in section 4. We consider a two-dimensional system of N soft spheres ofequal diameter σ . The soft sphere interactions are repulsive | ?? nd are modeled by a shifted and truncated (12,6) Lennard-Jones Potential U SS ( r ij ) = 4 (cid:15) (cid:0) ( σ/r ij ) − ( σ/r ij ) + 1 / (cid:1) (1)which is cut off at r c,SSij = 2 / σ . Here, r ij = | r j − r i | is theparticle center-to-center distance and (cid:15) sets the unit of energy.The crossed orthogonal external fields induce orthogonaldipole moments µ m = µ e m and µ e = µ e e which we termfor simplicity as ’magnetic’ and ’electric’ dipoles (althoughthe model is also appropriate for two electric moments). Thecoordinate frame is adjusted to coincide with the directionsof these moments so that e m = e x and e e = e y . In gen-eral these moments could have different absolute values butfor simplicity they are assumed to be equal. The two typesof dipole moments are also assumed to be independent fromeach other and interact only with dipole moments of the sametype on other particles. Intuitively, one would model the in-teraction energy between dipoles of particles 1 and 2 by thepoint-dipole potential U αdip ( r ) = µ α · µ α r − µ α · r )( µ α · r ) r , (2)where α indicates the dipole type as being either e or m . Dueto the constraint µ α (cid:107) µ α it follows that U αdip ( r ) = µ α µ α r (1 − r · e α ) r ) . (3)The resulting total dipolar interaction between two particles isthe sum of the dipolar potentials stemming from the magneticand electric dipoles, respectively. Using µ = | µ αi | and therelation ( r · e e + r · e m ) = r (which holds since µ e and µ m are orthogonal) we obtain U edip ( r ) + U mdip ( r ) = − µ r . (4)The resulting interaction on the right side of Eq. (4) is anisotropic, purely attractive interaction that lacks any kindof directional character. Therefore, the potential defined inEq. (4) can not generate any rectangular structures as observedin experiments . Moreover, having in mind that the field-induced spatial separation of charges in the ionic layer of asuspended particle is comparable to the particle diameter σ , apoint dipole model seems even more unreasonable.We therefore introduce an alternative model. Each dipolemoment µ α (with α = e, m ) is replaced by two oppositecharges − q α = q α which are shifted out of the particle cen-ter by a vector δ α k = ( − k δ e α , with k = 1 , . The vec-tor δ α k points either parallel ( k = 2 ) or antiparallel ( k = 1 )along the corresponding point dipole moment µ α . Indepen-dent of their shift or type, we set all charges to the same abso-lute value q = | q α k | = 2 . (cid:15)/σ ) − / . The choice of the value Fig. 2 (Color online) Normalized direction-dependent pair interac-tion U ( r ij ) [see Eq. (9)] between a particle in the center of the co-ordinate frame and a second particle (indicated as black circle in(a)) at various positions r ij for three different charge separations δ = 0 . , . , . σ corresponding to (a), (b), (c). Sterically ex-cluded areas are indicated by white circles. (d) Interaction energyat distance r ij = σ as function of φ , the angle measured in multiplesof π against the x-axis, for δ = 0 . σ (yellow), δ = 0 . σ (purple)and δ = 0 . σ (black). . is essentially arbitrary, as we will later normalize the in-teraction energy to eliminate the dependence of its magnitudeon the charge separation δ [see Eq. (8) below]. Charges k and l on different particles i and j interact via a Yukawa potential U α k α l ij ( r ij ) = − q exp ( − κr α k α l ij ) r α k α l ij (5)with r α k α l ij = | r j − r i + δ α l − δ α k | . The inverse screeninglength is choosen to κ = 4 . σ − and a radial cutoff r c = 4 . σ is applied. A schematic representation of the model with its in-ternal arrangement of ’charges’ is shown in Fig. 1. Due to theYukawa-like interaction our model lacks the long-range char-acter of true dipolar interactions. However, similar modelswith comparable interaction ranges have been used to describedipolar colloids in the context of discontinous molecular dy-namics simulations . We note that mimicking magneticdipoles via spatially separated ’charges’ appears somewhat ar-tificial from a physical point of view. However, here we couldthink of a particle with a strongly inhomogenous distributionof magnetic moments. Moreover, the idea behind our ansatz isnot to mimick a particular complex colloid, but rather to pro-vide a generic model for field-induced directional interactions.The arrangment of charges inside particles then results in a ?? | air-interaction U DIP ( r ij ) given by U DIP ( r ij ) = (cid:88) k,l =1 [ U e k e l ij ( r ij ) + U m k m l ij ( r ij )] . (6)In principle, U DIP ( r ij ) is a function of q and δ . To facilitatethe comparison between the interactions at different δ ( q ischosen to be constant), we normalize U DIP ( r ij ) according to ˜ U DIP ( r ij ) = U DIP ( r ij ) × u/U DIP ( σ e α ) (7)where the constant u = − . (cid:15) is calculated from the unnor-malized energy U DIP ( σ e α ) with model parameters δ = 0 . σ and q = 2 . (cid:15)/σ ) − / . This procedure ensures that the nor-malized energy between two particles at contact ( r ij = σ ) anddirection r ij = σ e α (pointing along one of the fields) has theconstant value u for all δ , that is ˜ U DIP ( σ e α ) = u. (8)The full pair interaction of our model is then given by U ( r ij ) = U SS ( r ij ) + ˜ U DIP ( r ij ) . (9)The resulting potential is illustrated in Fig. 2(a)-(c) for a parti-cle in the center of the coordinate frame and a second particleat various distances r ij and angles φ = arccos( r · e x /r ) with ’charge’ separations δ = 0 . , . , . σ . The value δ = 0 . σ is motivated by our simulation results presented inSec. 3.1. Sterically-excluded areas are shown in white and en-ergy values are color coded in units of (cid:15) . The weak anisotropyof the resulting particle interactions at small δ (where one es-sentially adds two dipolar potentials, see Eq. (4)) transformsto a patchy-like pattern by increasing δ . Energy minima be-come more and more locally restricted and particle interac-tions become directional in character. This is also seen inFig. 2(d) which gives the energy between two particles in con-tact as function of φ for different δ . From Fig. 2(d) we alsosee that, independent of the ’charge’ separation δ , the min-ima of the full interaction potential [see Eq. (9)] occur forconnection vectors r ij = σ e e and r ij = σ e m (i.e., pointingalong the fields). Note that this already holds for the unnor-malized energy given in. Eq. (6). Simulations are performedwith N = 1800 to particles at a range of reduced num-ber densities ρ ∗ = ρσ and temperatures T ∗ = k B T /(cid:15) , in asquare-shaped simulation cell with periodic boundary condi-tions. The equations of motion ˙ r i = T ∗ γ N (cid:88) j =1 ∇ U ( r ij ) + ζ i ( T ∗ ) (10)are solved via the Euler scheme with an integration stepwidth ∆ t = 10 − τ b , where τ b is the Brownian timescale defined by τ b = σ γ/T ∗ , γ = 1 . is a friction constant and ζ i ( T ∗ ) is aGaussian noise vector which acts on particle i and fulfills therelations (cid:104) ζ i (cid:105) = 0 and (cid:104) ζ i ( t ) ζ j ( t (cid:48) ) (cid:105) = 2 γk B T δ ij δ ( t − t (cid:48) ) .We perform simulations for up to τ b . To characterize the structure of the systems we consider sev-eral quantities. The first one is the mean coordination number ¯ z = 1 N N (cid:88) i =1 z i , (11)where z i is the number of neighbors of particle i and the sumis over all particles. In the following, two particles are consid-ered to be nearest neighbors if their center-to-center distanceis smaller than r b = 1 . σ .To identify local particle arrangements, the orientationalbond order parameter is of special importance. For particle k it is given by φ nk = 1 z k z k (cid:88) l =1 | exp( inθ klλ ) | (12)with z k being the number of neighbors and θ klλ = arccos( r kl · r kλ / ( r kl r kλ )) being the angle between the bond of particle k and its neighbor l measured against a randomly chosen bondof particle k to one of its neighbouring particles λ . Hence, φ nk = 0 for z k < . The integer value n determines the typeof order which is detected by this parameter. We concentrateon φ and φ to identify square (rectangular) and hexagonallattice types. Its ensemble average is calculated via Φ n = 1 N N (cid:88) i =1 φ ni . (13)The reversibility of ’bond’ formation and slipperyness of ex-isting bonds can be characterized by the bond and the bond-angle auto-correlation functions c b ( t ) and c a ( t ) . To evaluate c b ( t ) we assign a variable b ij ( t ) to each pair of particles ateach time step which is if the particles i and j are nearestneighbors or zero otherwise. The bond auto-correlation func-tion is then defined as c b ( t ) = (cid:104) b ij ( t ) b ij ( t ) (cid:105) , (14)where the brackets indicate an average over all pairs that arebonded at time t . The bond-angle auto-correlation function c a ( t ) is defined similarly by defining the unit vector a ij ( t ) = r ij ( t ) /r ij ( t ) , (15)such that c a ( t ) = (cid:104) − arccos ( a ij ( t ) · a ij ( t )) /π (cid:105) (16) | ?? ig. 3 (Color online) Simulation snapshots at ρ ∗ = 0 . and T ∗ = 0 . for (a) δ = 0 . σ , (b) δ = 0 . σ and (c) δ = 0 . σ . Particles arecolored according to their value of φ i . where we average again over all pairs. While c b gives theinformation on how stable bonds are over time, c a tells howstable their direction is over time. Note that in contrast to thetypical definition of correlation functions for stationary sys-tems , here the functions c b ( t ) and c a ( t ) are not independentof the time origin t .Finally, we consider the fractal dimension D f of parti-cle clusters, which is particularly important in the context ofDLCA. Clusters are defined as a set of particles with commonnext neighbors. The size of a cluster is then quantified by itsradius of gyration R g = 1 N cl N cl (cid:88) i =1 ( r i − ¯r ) , (17)where N cl is the number of particles in the cluster, and ¯r is thepostion of its center-of-mass. By plotting ln R g against ln N cl for different clusters, we extract the fractal dimension D f viathe relationship R g ∼ N /D f cl (see Ref. ). Our large-scale Brownian dynamics simulations show that thesystem is very sensitive to changes in temperature T ∗ , num-ber density ρ ∗ , and charge separation δ . In this large parameterspace we find a variety of different states ranging from smallfractal aggregates and single-chain structures at low tempera-tures to coarser, isolated or interconnected clusters at highertemperatures. In the following sections 3.1 - 3.3 we first dis-cuss the structure, the time correlation functions and the frac-tal dimensions at a low temperature and an intermediate den-sity, focussing on the impact of the model parameter δ . Insection 3.4 and 3.5 we then turn to the impact of temperatureand density. At first we study the system at low temperature T ∗ = 0 . andintermediate density ρ ∗ = 0 . for different charge separations δ . In Fig. 3 simulation snapshots for δ = 0 . σ, . σ, . σ at t = 300 τ b [see. Eq. (10) below] are shown, where τ b is theBrownian timescale. The colorcode reflects the orientationalbond order parameter φ i of each particle i . All three cases arecharacterized by clusters with irregular shapes. However, lo-cal particle arrangements differ strongly. While for δ = 0 . σ the particles aggregate in a hexagonal fashion, at δ = 0 . σ they aggregate into rectangular structures. At the intermedi-ate charge separation δ = 0 . σ , hexagonal order dominatesthe system; however, some clusters also reveal subsets of par-ticles in rectangular arrangements. A more quantitative de-scription is given by the orientational bond order parameters Φ shown in Fig. 4(a) as functions of δ . By increasing δ ,one observes a sharp transition at δ ≈ . σ from hexagonaltowards rectangular (square) order.The very presence of such a sharp transition can be ex-plained via energy arguments based on the δ -dependent pairpotential plotted in Figs. 2(a)-(d). To this end, we calculatethe energy U hexi ( δ ) = (cid:80) j =1 U ( r ij ) of a particle i with sixneighbors j , which are located in a hexagonal arrangement at’contact’ distance σ around i . Note that not all hexagonal con-figurations do have the same contact energy. This is due to theanisotropy of interactions, see Fig. 2(d). Therefore we con-sider a hexagonal configuration in which the contact energy isas low as possible (this configuration was found numerically).The dependence of this lowest contact energy U hexi ( δ ) on thecharge separation parameter is plotted in Fig. 5. Also shownis the corresponding energy U sqi ( δ ) = (cid:80) j =1 U ( r ij ) = 4 × u of a particle with four neighbors j located at distance σ in arectangular arrangement, i.e., in the energy minima around i ?? | ig. 4 (Color online) Results for simulations with N = 1800 attemperature T ∗ = 0 . and density ρ ∗ = 0 . . (a) Orientationalbond order parameters Φ for square (black) and Φ (yellow) forhexagonal particle arrangements. (b) Mean coordination number ¯ z as function of charge separation δ at times t = 100 , , τ b . (the quantity u was defined below Eq. (7)). Note that the en-ergy U sqi ( δ ) does not depend on δ according to Eq. (8). Asshown in Fig. 5, the two curves intersect at a ”critical” valueof δ = 0 . σ . Thus, the simple energy arguments alreadysuggest a transition between states with local hexagonal andsquare order, even though the predicted critical value is some-what larger than the value of δ = 0 . σ seen in the actualsimulations at finite temperature and density [see Fig. 4(a)].Further information is gained from the behavior of the meancoordination number as a function of δ plotted in Fig. 4(b) forthree different times t = 100 τ b , τ b and τ b . At all timesconsidered, ¯ z undergoes a steep decrease at δ ≈ . σ froma nearly constant value, ¯ z hex ≈ . , to a value ¯ z sq ≈ . .This behavior reflects, on the one hand, again the presenceof a sharp transition; on the other hand, the actual values of ¯ z hex (¯ z sq ) reveal the ”non-ideal” character of the aggregates interms of coordination numbers. For example, for δ > . σ we find that ¯ z and Φ decrease with δ , while Φ increases.However, this does not indicate a decline of the rectangu-lar order; it rather results from an increasing amount of par-ticles residing in chains oriented either in x- or y-direction.The coordination number z i of a particle i in such a chain is ≤ , leading to a mean coordination number ¯ z < . Further-more, the parameters φ i and φ i [see Eq. (12)] become unity Fig. 5 (Color online) Minimum energy of a particle with six neigh-bors in hexagonal arrangement as function of δ (black) and energyfor a particle in rectangular arrangement with 4 neighbors (red). for a particle forming exactly two bonds under an angle of π (straight chain). This does not affect Φ , which is alreadylarge at δ > . σ , but significantly increases Φ . Finally,the counter-intuitive decrease of Φ with δ results from theincreasing amount of particles with only one neighbor (e.g.,ends of chains appearing white in Fig. 3(c)). These particlesyield no contribution to Φ [see Eq. (12)].The ”non-ideal” values of ¯ z hex and ¯ z sq also explain whyour energy argument for the location of the hexagonal-to-square transition, which was based on ideal arrangements withsix and four neighbors, respectively, does yield the transitionvalue δ = 0 . σ rather than δ = 0 . σ obtained from sim-ulation. We can now reformulate the argument by using theactual mean coordination numbers extracted from our simu-lations, ¯ z sq = 3 . (instead of ) and ¯ z hex = 4 . (insteadof ). Following the calculations for the ideal arrangementsdescribed before, the energy of the square-like arrangement is U sq = 3 . × u . For the hexagonal arrangement, we use the av-erage minimum energy with either z i = 4 or z i = 5 neighbors,yielding ¯ U h (¯ z hex , δ ) = ( U hex (4 , δ ) + U hex (5 , δ )) / . The re-sulting critical value of the charge separation is δ ≈ . σ ,which coincides nicely with the transition value observed inour simulations. Although the local structures characterized by ¯ z and Φ per-sist, in general, over the simulation times considered, we arestill facing a transient (out-of-equilibrium) structure forma-tion as seen, e.g., from the slight increase of ¯ z with time inFig. 4(b). This raises a question about the typical ”lifetime” ofthe aggregates.To this end we now consider dynamical properties, namelythe bond and bond-angle auto-correlation functions, c b ( t ) and c a ( t ) . It is not reasonable to extract decay rates from thesefunctions (as it is usually done) because in transient states, de-cay rates are, strictly speaking, functions of time themselves.Still, it is interesting to see whether the temporal correlation | ?? ig. 6 (Color online) Time correlation functions obtained from sim-ulations with N = 1800 at temperature T ∗ = 0 . and density ρ ∗ = 0 . . (a) [(b)] Time evolution of the bond [angle] autocor-relation function c b ( t ) [ c a ( t ) ] for three different charge separations δ = 0 . σ, . σ and . σ colored in yellow, purple and black re-spectively. of bonds (bond angles) for different δ allows us to distinguishbetween qualitatively different aggregation regimes.Numerical results for c b ( t ) and c a ( t ) are plotted inFigs. 6(a) and (b), respectively, where we consider a large timerange up to t ≈ τ b . The time axis starts at the finite timewhen all the systems have formed stable aggregates. The datain Figs. 6(a) and (b) pertain to three representative values ofthe charge separation parameter related to the hexagonal struc-tures ( δ = 0 . σ ), rectangular structures ( δ = 0 . σ ), and to thetransition region ( δ = 0 . σ ). In the square regime ( δ = 0 . σ )the decay of both c b ( t ) and c a ( t ) is almost identical and veryslow. From this we conclude that the square regime is charac-terized by almost unbreakable bonds with fixed orientations.This is different in the hexagonal regime ( δ = 0 . σ ) where c b ( t ) remains nearly constant even after long times (meaningthat bond-breaking is very unlikely), while c a ( t ) decays muchfaster. Thus, the direction of bonds are less restricted. Weinterprete this behavior as evidence that two particles, thoughbeing bonded, are still able to rotate around each other to someextent. This is a characteristic feature of slippery bonds. Fi-nally, in the transition regime ( δ = 0 . σ ) both functions c b ( t ) and c a ( t ) decay significantly faster than in the other cases,with the decay of the bond-angle correlation function beingeven more pronounced. In that sense we may consider thebonds in the transition region also as slippery (although lesslong-lived than in the other cases).We conclude that the different structural regimes identifiedin the preceding section are indeed characterized by differentrelaxational dynamics. Moreover, all of the observed aggre-gates have lifetimes of at least several hundered τ b . Such long-lived bonds are indicative of diffusion limited cluster-clusteraggregation. In the next section we therefore consider the frac-tal dimension. Fig. 7 (Color online) Fractal dimension D f as a function of chargeseparation at ρ ∗ = 0 . and T ∗ = 0 . . At δ = 0 . σ we find abimodal distribution of fractal dimension with peaks at D f = 1 . (solid line) and . (dashed line). In Fig. 7 the fractal dimension D f is shown as a function of δ at time t = 250 τ b , density ρ ∗ = 0 . and temperature T ∗ =0 . . We find that D f increases slightly with δ but remainsin a range between . and . , except at δ = 0 . σ . There,the fractal dimension exhibits a bimodal distribution, takingvalues between D f ≈ . and D f ≈ . (dashed line inFig. 7).Despite these variations, the values of D f found here aresignificantly smaller than the fractal dimension D f = 1 . ob-served in earlier studies of DLCA in two-dimensional contin-uous (off-lattice) systems . Except for the case δ = 0 . σ ,the values in Fig. 7 are comparable with previous findings forDLCA in two-dimensional lattice systems and systems withspatial or interaction anisotropies . The present system isindeed anisotropic in the sense that the external fields imposepreferences on the directions of particle bonds and thereforealso on the orientations of aggregates. This effect is mostpronounced in the rectangular regime ( δ = 0 . σ ). There-fore, it is plausible that our system undergoes a special caseof anisotropic DLCA, in (quantitative) accordance with ex-perimental results and theoretical predictions . Weshould note that, due to our simulation method, the clustersizes (typically involving − particles) are relativelysmall compared to the particle numbers considered in the lit-erature ( particles) . A more detailed analysis of theimpact of anisotropic interactions on the fractal structure isbeyond the scope of this study.We also relate our findings to the newer concept of slip-pery DLCA , where the bonds are essentially unbreakablebut able to rotate. Indeed, as discussed in section 3.2, bondsare slippery in nature for small δ in the hexagonal regime.For three-dimensional systems it has been reported thatthe fractal dimension D f remains the same for slippery andclassical DLCA, while the mean coordination number ¯ z dif-fers. Specifically, ¯ z is significantly higher for sDLCA . ?? | ig. 8 (Color online) Temperature dependence of the system properties at density ρ ∗ = 0 . for charge separations δ = 0 . , . , . σ coloredin yellow, purple and black, respectively. (a) Fractal dimension D f evaluated at t ≈ τ b , (b) Mean coordination number, (c) Orientationalorder parameters Φ and Φ . The same observation emerges when we consider our valuesof ¯ z plotted in Fig. 4(b), from which one sees a pronounceddecrease of ¯ z upon entering the square (DLCA) regime. How-ever, in contrast to earlier studies we find D f to slightly in-crease with δ , especially in the hexagonal regime. We in-terpret this behavior as a consequence of the fact that bind-ing energies in the hexagonal regime decrease with increasingvalues of δ , while they remain constant in the square regime(see Fig. 5). The corresponding stability of bonds should becorrelated to the binding energies which explains the slightlyincreasing values of D f in the hexagonal regime.Finally, in the transition region ( δ = 0 . σ ) we found a bi-modal distribution of the fractal dimensions D f with maximaat D f ≈ . and D f ≈ . . This second maximum corre-sponds to only ≈ of the considered cases (twelve inde-pendent simulation runs). The first maximum at D f ≈ . therefore clearly dominates and fits nicely to the functionaldependence of D f on δ (see Fig. 7). We assume that theless frequent peak results from a switching of the local struc-tures between hexagonal and rectangular arrangements, whichis accompanied by a significantly larger bond-breaking prob-ability (see Fig. 6(c)). Again this allows compactification ofaggregates and increases the fractal dimension in the transitionregime. Diffusion limited aggregation is restricted to systems with at-tractive particle interactions much stronger than k B T . By in-creasing the temperature sufficiently, thermal fluctuations be-come able to break bonds which results in a faster decay of thethe bond auto-correlation functions and a compactification ofaggregates. Indeed, for square lattice models it was found that D f is a monotonically increasing function of temperature .In Fig. 8(a) the fractal dimension D f of the present model isplotted as a function of temperature T ∗ for charge separations δ = 0 . σ, . σ and . σ . We first concentrate on the case δ = 0 . σ , correspondingto the square regime at low T ∗ . In the range of very low tem-peratures T ∗ < . , the fractal dimension is small and staysessentially constant. Increasing T ∗ towards slightly larger val-ues then leads to an increase of D f , reflecting the (expected)compactification. This increase of D f is accompanied by anincrease of the mean coordination number ¯ z [see Fig. 8(b)]within the temperature range considered, indicating the grow-ing number of bonds due to local and global structural re-configurations. The corresponding changes in the stability ofthe bonds are illustrated in Fig. 9, where we have plotted thetime evolution of c b ( t ) for several temperatures (at δ = 0 . σ ).Clearly, the decay of c b ( t ) becomes faster for higher tempera-tures. This is the reason why structural reconfigurations and,in consequence, compactification of aggregates becomes pos-sible.These trends persist until T ∗ f,sq ≈ . , beyond which thesystem at δ = 0 . σ starts to behave in a qualitatively differ-ent way. The mean coodination number ¯ z displays a maxi-mum and subsequently a rapid decay. We also find that thefractal dimension has not yet reached its maximum value at T ∗ f,sq ; this maximum occurs at the slightly larger temperature T ∗ ≈ . (see Fig. 8(a)). This ’delay’ of D f can be under- Fig. 9 (Color online) Bond auto correlation function c b ( t ) for dif-ferent temperatures T ∗ at charge separation δ = 0 . σ and density ρ ∗ = 0 . . | ?? ig. 10 (Color online) Simulation snapshots at ρ ∗ = 0 . with δ = 0 . σ at (a) T ∗ = 0 . , (b) T ∗ = 0 . , and (c) T ∗ = 0 . . Particles arecolored according to their value φ i . stood from the fact that, upon the entrance of bond-breaking,filigree parts of the aggregates are more likely affected thanmore compact ones. Hence, the fraction of ’compact’ smallaggregates still grows. Even more important, the function Φ ( T ∗ ) in Fig. 8(c) displays a pronounced decay of rectan-gular order for T ∗ > T f,sq . From the sum of these indicationswe conclude that, at temperatures higher than T ∗ f,sq ≈ . ,the system transforms into a (stable or metastable) fluid phase.In this fluid phase, the overall structure starts to become ho-mogeneous and isotropic, while the local structures involveonly a small number of bonds with short bond-life times.For the system at δ = 0 . σ (hexagonal structure at low T ∗ ),an estimate of the ”fluidization” temperature T ∗ f,hex based onthe behavior of order parameters, coordination number andfractal dimension is more speculative. Nevertheless, the datasuggest that T ∗ f,hex > T ∗ f,sq . This is indicated, first, by thefact that Φ ( T ∗ ) decays only very slowly with temperatureuntil T ∗ ≈ . (see Fig. 8(c)). Second, the mean coordinationnumber shows only a weak maximum (and no fast decay after-wards) compared to the case δ = 0 . σ . Third, the fractal di-mension keeps increasing with T ∗ for all considered tempera-tures T ∗ < . . Therefore we conclude that T ∗ f,hex > . . Weunderstand this higher fluidization temperature at δ = 0 . σ from the fact that binding energies in hexagonal structures arelarger; therefore, higher coupling energies must be overcome.To further justify these interpretations, particularly theemergence of fluid phases, we performed a stability analy-sis of the homogenous isotropic high temperature state basedon mean-field density functional theory (DFT). Specifically,we consider the isothermal compressibility χ T . Positivevalues of χ T imply that the homogeneous (fluid) phase isstable, whereas negative values indicate that this phase isunstable. Specifically, the instability arises against long-wavelength density fluctuations, i.e. condensation. According to Kirkwood-Buff theory one has χ − T ∝ − ρ ˜ c ( k = 0) , (18)where ˜ c (0) is the Fourier transform of the direct correlationfunction (DCF) c ( r ) in the limit of long-wavelengths ( k → ). We approximate the DCF for distances r ij > σ accordingto a mean field (MF) approximation, that is c MF ( r ) = − ( k B T ) − U ( r ) , r > σ, (19)and use the Percus-Yevick DCF c HS ( r ) of a pure hard-sphere fluid for | r | ≤ σ . The full DCF is then given by c ( r ) = c HS ( r ) + c MF ( r ) . (20)In Fig. 11 we present numerical results for the expression − ρ ˜ c (0) at ρ ∗ = 0 . as function of temperature. At low T ∗ ,all systems are characterized by negative values of − ρ ˜ c (0) .This indicates that the homogeneous isotropic phase is unsta-ble, consistent with the results of our simulations. Upon in-creasing T ∗ the mean-field compressibility χ T then becomes Fig. 11 (Color online) Numerical solutions to Eq.18 as function of T ∗ for density ρ ∗ = 0 . and charge separations δ = 0 . , . , . σ colored in yellow, purple and black, respectively. ?? | ndeed positive for all charge separations considered. Specif-ically, for δ = 0 . σ the change of sign (related to a ”spinodalpoint”) occurs at T ∗ f,sq = 0 . and for δ = 0 . σ at the muchhigher temperature T ∗ f,hex = 0 . . These values are in surpris-ingly good agreement with our estimates for the ”fluidization”temperatures based on the order parameter plots.The case δ = 0 . σ is again different. Here we find[see Fig. 8(b)] that, starting from low temperatures inside theDLCA regime, the mean coordination number monotonicallydecreases. However, this does not indicate ”fluidization” butrather a gradual transition from a state with dominant hexago-nal order towards a mixed state comprised of coexisting clus-ters with local hexagonal and square-like order. Indeed, [seeFig. 8(c)], the orientational order parameters Φ and Φ revealthat the fraction of particles bound in square clusters increaseswith T ∗ and finally overtakes the fraction of particles involvedin hexagonal clusters at T ∗ ≈ . . Corresponding snapshotsof simulation results are shown in Fig. 10. At all temperaturesconsidered one observes separated clusters. With increasingtemperature their shape becomes more regular, while the lo-cal rectangular order becomes more pronounced. Finally, at T ∗ = 0 . the fractal dimension D f and the square orderparameter Φ reach their maximum values, suggesting a ”flu-idization” similar to the behavior observed at other values of δ .Interestingly, our stability analysis [see Eq. (18)] indicates aninstability at the same temperature T ∗ f = 0 . . With this sur-prisingly accurate agreement between theory and simulation,we conclude that in the transition regime ( δ = 0 . σ ), increas-ing thermal fluctuations first push the system from a domi-nantly hexagonal state into a rectangular one, which then en-ters a metastable fluid phase after passing the ”spinodal point”. In this section we revisit the system behavior at the low tem-perature T ∗ = 0 . , but consider different densities in therange ρ ∗ ≤ . . Whereas low-density systems at T ∗ = 0 . display DLCA as discussed in sections 3 A-C, this aggrega-tion mechanism is expected to disappear at higher densities:here, the particles are just unable to diffuse sufficiently freely.Rather, the particles will very frequently collide and then im-mediately form rigid bonds. A typical structure at the high-est density considered, ρ ∗ = 0 . , and separation parameter δ = 0 . σ is shown in Fig. 12. Clearly, the system is per-colated, that is, the particles form a single, system-spanningcluster. Interestingly, this cluster is composed of extended re-gions characterized by either square-like order or hexagonalorder. We note that, at δ = 0 . σ , simultaneous appearanceof clusters with both types of order also occurs at low den-sities and higher temperatures (see section 3.4). However, atthe high density considered here the regions of each type arelarger and the particle arrangements are much more regular Fig. 12 (Color online) System at T ∗ = 0 . , density ρ ∗ = 0 . and δ = 0 . σ . The colorcode gives the orientational bond-orderparameter φ i of each particle i . (i.e., there are less defects).To better understand the impact of the density on the clusterstructures we plot in Fig. 13(a) the orientational bond orderparameters Φ and Φ as functions of ρ ∗ for δ = 0 . σ (at δ = 0 . σ and δ = 0 . σ the order parameters are essentiallyindependent of the density). From Fig. 13(a) it is seen thatthe amount of rectangular (hexagonal) order sharply increases(decreases) at a density of ρ ∗ ≈ . . This is a surprisingresult as one would expect that, upon compressing the system,close-packed, hexagonal structures rather become more likely.However, at the low temperature considered here, structuralreorganization is strongly hindered.We also note that all of the systems investigated at densi-ties ρ ∗ > . turned out to be percolated (suggesting that thevalue ρ ∗ = 0 . is indeed related to the percolation transi-tion). It thus seems that the percolation tends to stabilize theinitially formed square-lattice symmetry, as the subsequent re-organization is hindered by the lack of mobility. In effect, weare faced with quenched states that could not densify withinthe time domain studied. This interpretation is also consistentwith the decrease of the mean coordination number once thesystem is percolated ( ρ ∗ > . ) as shown in Fig. 13(b). In this work we propose a new model for field-directed ag-gregation of colloidal particles with anisotropic interactionsinduced by external fields. The model was inspired by recentexperimental work on novel colloidal particles inwhich external fields can induce two, essentially decoupled,dipoles. |1– ?? ig. 13 (Color online) (a) Orientational bond order parameter Φ and Φ as function of density ρ ∗ for δ = 0 . σ at T ∗ = 0 . . (b)Mean coordination number ¯ z as function of ρ ∗ at T ∗ = 0 . fordifferent δ = 0 . σ, . σ and . σ colored in yellow, purple andblack, respectively. The formation of particle networks with multiple percola-tion directions can find application in a range of new materi-als with anisotropic electrical and thermal conduction, mag-netic or electric polarizability or unusual rheological proper-ties. The aggregated clusters can be dispersed in liquid, whilethe percolated networks can be embedded in a polymer or gelmedium . The key to the fabrication of such novel classesof materials containing particle clusters and networks is thecontrol of the process parameters to obtain the desired inter-connectivity, density and structure.Against this background, the focus of our theoretical studywas to understand the formation of transient , aggregated struc-tures appearing at low-temperatures. Performing large-scaleBD simulations we have found that, depending on the distri-bution of the field-induced attractive ”sites” in the particles,different aggregation mechanisms arise. These have been an-alyzed via appropriate structural order parameters, bond time-correlation functions as well as by the fractal dimension. OurBD results demonstrate that by varying the charge separationparameter, that is, the distribution of attractive sites, the sys-tems transform from DLCA (essentially rigid bonds) towardssDLCA (slippery bonds). Moreover, we show that the changeof aggregation behavior is accompanied by significant changesof the local cluster structure.Indeed, the cluster structure can be easily manipulated byexploiting the interplay between temperature, density andmodel parameter δ . This allows formation of unexpectedstructures e.g., pronounced rectangular packing instead ofclosed packed hexagonal structures by increasing density.This unusual behavior appears to be dictated by the inabil-ity of the originally formed lattices with square symmetryto re-arrange into more dense hexagonal lattices. It has po-tentially important consequences for colloidal assembly, as ispoints out the ability to use mutidirectional field-driven assem-bly for the making of lower-density, yet highly interconnected,phases. Future research should focus on a more detailed investiga-tion of the interplay between the aggregation mechanisms ob-served here (anisotropic and slippery DLCA), and the equi-librium phase behavior, particularly the location of a conden-sation transition and of percolation at higher densities. Thisincludes investigation of the influence of entropy which wedid not discuss but is expected to strongly influence the aggre-gation behavior .In conclusion, our study makes an important contribution tocurrent research on structural and dynamical phenomena ac-companying self-assembly of complex colloids . In partic-ular, we have introduced a generic model describing colloidsin multidirectional fields yielding tunable multipolar interac-tions. Our study thus provides a microscopic understanding ofaggregation processes in such systems. Acknowledgements
We gratefully acknowledge stimulating discussions with Bhu-vnesh Bharti (NCSU). This work was supported by the NSF’sResearch Triangle MRSEC, DMR-1121107 and by The Ger-man Research Foundation via the International ResearchTraining Group (IRTG) 1524 ’Self-Assembled Soft MatterNano-Structures at Interfaces’.
References
Soft Matter , 2008, , 663.2 Yufeng Wang, Yu Wang, Dana R. Breed, Vinothan N.Manoharan, Lang Feng, Andrew D. Hollingsworth, MarcusWeck1, and David J. Pine, Nature , 2012, , 513 M. J. Salomon,
Current Opinion in Colloid & Interface Sci-ence , 2011, , 158.4 S. Scanna, L. Rossi and D. J. Pine, J. Am. Chem. Soc. , 2012, , 6112.5 F. Ma, D. T. Wu and N. Wu,
J. Am. Chem. Soc. , 2013, ,7839.6 J. Yan, M. Bloom, S. C. Bae, E. Luijten and S. Granick,
Nature , 2012, , 578.7 G. Rosenthal, K. E. Gubbins, and S. H. L. Klapp,
J. Chem.Phys. , 2012, , 174901.8 S. O. Lumsdon, E. W. Kaler, and O. D. Velev
Langmuir ,2004, , 2108.9 B. Bharti and O. D. Velev, Langmuir , 2015.10 A. Ruditskiy, B. Ren and I. Kretzschmar,
Soft Matter Langmuir , 2008, , 13312.12 S. Gangwal, A. Pawar, I. Kretzschmar and O. D. Velev, Soft Matter , 2010, , 1413. ?? | Phys. Rev. E , 2014, ,013105.14 H. L¨owen, J. Phys.: Condens. Matter , 2008, , 404201.15 S. Elborai, D.-K. Kim, X. He, S.-H. Lee, S. Rhodes, andM. Zahn, Journal of Applied Physics , 2005, , 10Q303.16 A. Dawar and A. Chandra, Physics. Letters. A. , 2014, ,2951.17 H. Schmidle, S. J¨ager, C. K. Hall, O. D. Velev and S. H.L. Klapp,
Soft Matter , 2013, , 2518.18 H. Schmidle, C. K. Hall, O. D. Velev, and S. H. L. Klapp, Soft Matter , 2012, , 1521.19 J. Russo, J. M. Tavares, P. I. C. Teixeira, M. M. Telo daGama and F. Sciortino, Phys. Rev. Lett. , 2011, , 085703.20 E. Zaccarelli,
J. Phys.: Condens. Matter. , 2007, ,323101.21 S. Correzzi, C. D. Michele, E. Zaccarelli, P. Tartaglia, andF. Sciortino, J. Phys. Chem. B. , 2008, , 1233.22 T. A. Witten Jr. and L. M. Sander,
Phys. Rev. Lett , 1981, , 1400.23 P. Meakin, Physica D , 1995, , 104.24 P. Meakin, Phys. Rev. A , 1988, , 4799.25 S. Tolman and P. Meakin, Phys. Rev. A. , 1989, , 428.26 R. Zhang, P. K. Jha, and M. O. Cruz, Soft Matter , 2013, ,5042.27 M. D. Haw. M. Sievwright, W. C. K. Poon, and P. N.Pusey, Advances in Colloid and Interface Science , 1995, ,1.28 S. Babu, J.C. Gimel, and T. Nicolai, Eur. Phys. J. E. , 2008, , 297. 29 C. R. Seager, and T. G. Mason, Phys. Rev. E. , 2007, ,011406.30 G. Helgesen, A. T. Skjeltrop, P. M. Mors, R. Botet, and R.Julien, Phys. Rev. Lett. , 1988, , 1736.31 N. Yoshioka, I, Varga, F. Kun, S. Yukawa, and N. Ito, Phys.Rev. E. , 2005, , 061403.32 J. Byrom, and L. Biswal, Soft Matter , 2013, , 9167.33 M. Suzuki, F. Kun, and N. Ito, Phys. Rev. E , 2009, ,021402.34 M. N. Popescu, H. G. E. Hentschel, and F. Family, Phys.Rev. E. , 2004, , 061403.35 R. C. Ball, Physica , 62.36 J.-M. Jin, K. Parbhakar, L. H. Dao, and K. H. Lee,
Phys.Rev. E. , 1996, , 997.37 B. Bharti and O. D. Velev, to be published .38 F. L. Bragaa, and M.S. Ribeirob, Computer Physics Com-munications , 2011, , 1602.39 A. Goyal, C. K. Hall, and O. D. Velev,
Phys. Rev. E , 2008, , 031401.40 D. L. Ermak, J. Chem. Phys. , 1975, , 4189.41 J.-P. Hanson and I. R. McDonald, Academic Press , 2006.42 J. G. Kirkwood and F. P. Buff,
J. Chem. Phys. , 1951, ,774.43 X. Guo and U. Riebel, J. Chem. Phys. , 2006, , 144504.44 P. J. Krommenhoek and J. B. Tracy,
Part. Part. Syst. Char-act. , 2013, , 759.45 X. Mao, Q. Chen, and S. Granick, Nature Materials , 2013, , 217. |1–