Self-assembly and clustering of magnetic peapod-like rods with tunable directional interaction

Based on extensive Langevin Dynamics simulations we investigate the structural properties of a two-dimensional ensemble of magnetic rods with a peapod-like morphology, i.e, rods consisting of aligned single dipolar beads. Self-assembled configurations are studied for different directions of the dipole with respect to the rod axis. We found that with increasing misalignment of the dipole from the rod axis, the smaller the packing fraction at which the percolation transition is found. For the same density, the system exhibits different aggregation states for different misalignment. We also study the stability of the percolated structures with respect to temperature, which is found to be affected by the microstructure of the assembly of rods.


Introduction
Many efforts are currently devoted to the search of new functionalized particles in order to satisfy the constant need for materials with different properties. Recent advances in particle synthesis has resulted in a rapid growth of this field of material science. Colloids with directional interactions are promising candidates [1][2][3]. By tuning the direction, the self-assembling process can be controlled leading to the emergence of specific structures. A paradigm example are dipolar colloids whose interactions are governed by permanent or field-induced, magnetic or electric dipole moments, as well as particles with more complex multipolar interactions. Nanoparticles (NP) with a magnetic mono-domain (MN), are a typical example with a wide range of applications, including magnetic fluids [4], biomedicine [5], Magnetic Resonance Imaging (MRI) [6], data storage [7], among others. Currently growing attention is being paid to magnetic particles with anisotropic shape, due to their more complex properties when compared to those with spherical shape, as for example magnetic birefringence [8] and thermal conductivity [9]. Beyond anisotropy in the shape of the particles, the same attention was also addressed to cases where the anisotropy is in the location of the dipole with respect to the center of symmetry of the particle. In recent theoretical works the structure of fluids containing spherical particles with embedded off-centered magnetic dipoles [10,11] was investigated. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Beyond the nature of the interaction itself, the morphology of the NP can be used as a controlling parameter to functionalize the MN through the interaction direction, i.e. by tuning the dipole moment's direction of a peapod-like rigid rod. This is experimentally realised through the synthesis of monodisperse magnetoresponsive rods of desired diameter, length, and magnetic susceptibility [12].
Driven by the interest in the phase behavior of polar liquid crystals, earlier models of rodlike particles with a single longitudinal (or transversal) dipole moment have been intensively studied both by theory and numerical simulations [13][14][15][16]. We present here a numerical study of the self-assembly of a two dimensional system of stiff peapod-like straight filaments, composed of single magnetic dipolar beads that are rigidly linked one by one. A similar system was studied earlier by experiment [17] and simulations [18,19]. The interaction is given by the superposition of the dipolar fields of each dipole bead. Differently of our previous work [19], in this study the direction of the dipole moment is altered with respect to the rod axis, opening the possibility of a plethora of new kinds of assembled clusters. Many experiments involving assemblies of colloids are actually done at surfaces and/or thin films [20][21][22][23][24][25], which motivates us to explore two dimensional (2D) systems.
We also study the connectivity properties of the present system by focusing on the percolation transition as a function of the direction of the dipole moment. The percolation transition is related to gelation in attractive-driven colloidal systems. Percolation behavior is also of great relevance in highly connected materials due to the possibility of enhancing the electrical and thermal conductivity [17,26]. In addition, there is a relation with the change of the viscosity in systems with sufficient strong bond strength [27,28].
The interaction strength may either increase or decrease the volume fraction required for percolation according to the definition of the separation at which two particles are considered to be connected, the dimensionality, and the proximity to the critical temperature [29]. As observed by Miller and Frenkel [30] for a system with adhesive hard spheres (AHS), the localization of the percolation threshold, the critical value of a given parameter where the percolation statistically happens, also depends on the interaction strength between particles, which is directly related with the temperature of the system. Therefore it is worth to explore possible routes under which the percolation transition is enhanced. Concerning the dimensionality, and, by using elongated particles, it was shown, in the 3D case, that an increase of the aspect ratio decreases the percolation threshold [18,31], while the opposite behavior is observed in the 2D case [19]. In the present work we study the equilibrium configurations and the percolation transition as a function of the angular misalignment of the dipole moment with respect to the rod axis, and we show that larger values of this misalignment yields an enhancement of the percolation transition, i.e., an infinite extended cluster is formed for lower density.
The paper is organized as follows: our model system and simulation details are given in Sect. 2. The obtained different cluster configurations are presented in Sect. 3. The connectivity properties are discussed in Sect. 4. Our conclusions are given in Sect. 5.

Model and simulation methods
Extensive Langevin Dynamics simulations were performed to study a two-dimensional (2D) system consisting of typically, unless stated otherwise, N = 840 identical stiff rods of aspect ratio l = 3. The phase behavior of a mono-dispersed system with the same aspect ratio was recently studied [32] and is considered an established reference system [33]. For suspensions studied experimentally, the aspect ratio l = 3 is in the lowest accessible limit [34]. The magnetic nature of the rod is simulated by attaching a point dipole of permanent magnetic moment μ at the center of each bead (see Fig 1). The orientation of the dipoles with respect to the axial direction of the rod is given by the angle C, as illustrated in Fig 1. To model the dipolar particles we use a dipolar soft sphere (DSS) potential [19], consisting of the repulsive part of the Lennard-Jones (LJ) potential u rep and a point-like dipole-dipole interaction part u D . The total interaction energy between rods a and b is the sum of the pair interaction terms between their respective dipolar spheres (DS): where: with σ the diameter of each bead, and is the LJ soft-repulsion constant, R a,b = R b − R a is the vector which connects the center of rod b with the center of rod a. The orientation of rods a and b are given by θ a and θ b , respectively. The vector r a;b jm connects the center of bead m of rod b with the center of bead j of rod a (see Fig 1). The force on bead m due to bead j is given by: The torque on bead m is [19]: where d m is the vector connecting the center of bead m (rod b) with the center of rod b as illustrated in Fig 1, and B jm is the magnetic field generated by the dipole moment μ j at the position of the dipole μ m , which is given by: The summations in Eq (6) are considered only for dipoles belonging to distinct rods. The orientation of the rods is given by the unitary vector s given by s = d m / j d m j. The translational and rotational Langevin equations of motion of rod b with mass M b and moment of inertia I b , are given by: where is the angular velocity, F b and N b are the total force and torque acting on rod b, respectively, while Γ T and Γ R are the translational and rotational friction constants. ξ T b and ξ R b are the Gaussian random force and torque, respectively, which obey the following white noise conditions: For rodlike particles, the translational friction constant is a combination of the parallel and perpendicular components with respect to the rod axis, so that the total translational diffusion coefficient [35]. However, since the dynamical properties are not studied in the present work, for equilibrium simulations the values of Γ T as well as Γ R are irrelevant. We introduce the ratio Γ T / Γ R = 4/3 because such a ratio is already known to produce fast relaxation to equilibrium in similar systems [35]. We define the reduced unit of time as t Ã ¼ t= , where M is the mass of the rod. The energy is given in reduced units as U Ã = U/, the dipole moment in dimensionless units as , and the dimensionless distances as of r Ã = r/σ. Unless stated otherwise, the ratio of the thermal energy to the soft-sphere repulsion constant is chosen to be k B T/ = 0.1, where /k B is the temperature unit and k B is the Boltzmann constant. Periodic boundary conditions are taken in both spatial directions. Since the dipolar pair interaction falls off as (r −3 ), we take the simulation box sufficiently large such that no special long-range summation techniques [36] are needed. We define the packing fraction as η = N beads π(σ/2) 2 /L 2 , where N beads = 2520 is the total number of dipolar beads of the system and L 2 is the simulation box area. Since N beads = lN, we can rewrite the packing fraction as η = ρ Ã lπ/4, where ρ Ã is the dimensionless density ρ Ã = ρσ 2 , and ρ = N/L 2 . The reduced time step is typically in the range δt Ã = 10 −4 − 10 −3 .
All simulations are performed at a very slow cooling rate. Initially, we set the temperature at k B T/ = 2 and slowly reduced it in steps Δk B T/ = 0.05, each 5 × 10 5 time steps, till the final temperature is reached. The quantities of interest are then averaged over more than 10 6 time steps. All the beads from all rods have the same dipole moment whose magnitude we set as μ Ã = 1. Common experimental values of μ Ã2 at room temperature ranges in the interval 0.1 μ Ã2 10. For example, in experiments [37] carried out using aqueous dispersions of superparamagnetic microspheres of ferrite grains (Estapor (R) from Merck-reference M1-030/40) for r % 205 nm and M s % 6 × 10 4 A/m, the magnetization (M) of the particles is completely reversible and adjustable by an external magnetic field. If we consider T = 293 K and M % 22, 6% of M s , we obtain μ Ã % 1.

Cluster formation
We first present the dependence of the DSS pair interaction potential on the separation between rods. The study of the pair interaction potential is needed in order to understand the nature of the resulting many-body interaction and to help us to set the values of the useful parameters that help to understand the self-assembled structures. The dependence of the DSS pair interaction potential as a function of the separation between rods and minimized with respect to some characteristic angles is presented in Fig 2. Rather different potential profiles are obtained with changing values of C. For low values of C ( 30˚) the minima are located at r 0 /σ % 3, which corresponds to the aspect ratio of the rods. The values of α and θ which minimize the pair-interaction energy indicate that rods are favourably in the head-to-tail bond. As C is increased the position of the global minimum is displaced to smaller values of r 0 /σ, suggesting that the head-to-tail bond disappears, giving rise to the ribbon-like bond configuration (see Fig 1). Ribbon-like configurations are defined here as a side-by-side assembly as a consequence of the head-to-tail tendency of alignment between dipoles of beads in different rods and which are sufficiently displaced from the axial direction (C > 45˚). Similar ribbon-like configurations were very recently observed for microscopic magnetic ellipsoids [38,39], and in peanut-shaped colloids [40]. For C = 45˚we have a mixing of head-to-tail and ribbon-like arrangements with intermediate bonding energy and separation.
In order to find the separation (δ c ) used as a definition when two rods are bonded, we analyze the interrod separation related to the minimal energy value. From Fig 2, we observe that the largest interrod separation related to the global minimum is located around % 3.4σ for C = 15˚and 1.4σ for C = 90˚. In the former the rods are in the head-to-tail arrangement, while in the latter they are in the ribbon-like configuration. In both cases, the shortest separation between beads of different rods is % 1.4σ (bead-to-bead center distance). Therefore, we define here that two rods are bonded if the shortest separation between them is 1.4σ.
Since the attraction between magnetic rods becomes stronger for larger C (Fig 2), we expect that the ribbon-like configurations will become more stable, implying that the formation of clusters is facilitated in the many-body case. To show that this is indeed the case, we analyze the degree of polymerization [41], defined as: where N c is the number of clustered rods and N is the total number of rods. The polymerization as a function of the packing fraction η for different C is presented in Fig 3. There is a clear distinction of η-dependence of the polymerization for distinct values of C. For C 30˚the polymerization increases with increasing η, while for C > 30˚, the value of C does not change as η is increased. In principle we would expect that the tendency for aggregation should be stronger as η is increased. However, such a behavior is found only when C 30˚, where the head-to-tail arrangement is mostly observed. For C = 45˚, the configurations are a mixture of head-to-tail and ribbon-like arrangements (minimum energy at r/σ % 2, Fig  2) and the polymerization does not change significantly as η is changed. On the other hand, the larger C, the larger is F, since, based on the minimized pair interaction function (Fig 2), there is an increase of interrod attraction specially when the rods form a ribbon-like arrangement. For C ! 60˚, the value of F does not change as C is increased, and all the rods are connected to each other mostly arranged in the ribbon-like arrangement as the minimum energy is found for r/σ % 1.4. Self-assembly of magnetic rods with tunable directional interaction Some representative equilibrium configurations are presented in Fig 4 for packing fraction η = 0.1 and η = 0.3, and distinct C. We observe that the head-to-tail arrangements are found for C = 15˚and C = 30˚, while the ribbon-like configurations are observed for C = 60˚and C = 90˚. Such arrangements of the rods can be better characterised by computing the pair correlation function [42]: where R ab is the separation between the center of the rods a and b (see Fig 1). We present in Fig 5 the pair correlation function for η = 0.1 and for different values of C. As shown, the position of the peaks changes towards smaller r for larger C. For C = 15˚and C = 30˚the first largest peak is at r/σ % 3, which coincides with the value of the aspect ratio of the rods and is associated with the head-to-tail alignment. When C increases, multiple peaks proportional to r/σ 3 are found, and this is a consequence of the tendency to form side-byside arrangements, e.g., for C ! 60˚the first largest peak is at r/σ % 1.4. For intermediate values of the first largest peak position, i.e. r/σ % 2 for C = 45˚, is a consequence of the mixed structures of head-to-tail and ribbon-like alignments, as expected by the pair interaction in Fig  2. For C 45˚the g(r) does not maintain a constant structure and it loses the long-range correlation, which is a typical behavior of liquid-like structures. On the other hand, for C ! 60t he g(r) has regular peaks and long-range correlation, which is a typical behavior of solids. These results indicate, qualitatively, that the system goes from a liquid-like configuration (C 45˚) to a solid-like ribbon-like configuration (C ! 60˚), as a consequence of the fact that the attraction between rods in the ribbon-like configuration is much stronger than that observed in the linear head-to-tail arrangement, see Fig 2. In addition, such ribbon-like arrangements are expected to be more stable against thermal fluctuations.

Connectivity properties
Now we discuss the connectivity properties of the self-assembled structures. We focus our analysis on the study of the formation of percolated aggregates, which are defined as infinite connected clusters spanned over the system. The percolation transition is defined in the thermodynamic limit, where the average cluster size diverges [43]. We say that a configuration is percolated when, accounting for periodic boundary conditions, there is a percolating path [44], i.e., a cluster connected through opposite borders of the simulation box. A common feature of systems consisting of interacting particles subject to thermal fluctuations is that their bonds are transient. For low temperature, the lifetime of the bonds is sufficiently long, and the clusters are well defined over time. In our simulations, the temperature is one order of magnitude smaller than the average pair-interaction energy between rods, and the clusters are rather stable. In order to characterize the formation of the percolated clusters, we calculate an order parameter [45,46] defined as the fraction of monomers in the largest clusters, S max , i.e., where N larg is the number of rods belonging to the largest cluster. Previous studies showed that for finite size systems, the percolation transition is well characterized when the size of the largest cluster is at least 50% of the total number of particles, i.e., S max = 0.5 [45,46]. Such an order parameter is useful to study the connectivity properties of the system, specially gelation [45,47]. Here, we evaluate S max as a function of the packing fraction and the orientation of the dipole moments with respect to the axial direction of the rod. S max as a function of the packing fraction is presented in Fig 6(a) for different C. The percolation transition is shifted towards smaller packing fraction with increasing C, due to the stronger attraction between rods observed in these cases. As a consequence, the emergence of ribbon-like structures is facilitated. This is shown more clearly in Fig 6(b), where S max is presented as a function of C for different packing fractions. In all cases S max increases with increasing C, indicating that the rods interact more attractively, facilitating the formation of larger clusters. For small enough η (≲ 0.1) there is no formation of an infinite cluster independent of C. Also, we observe that S max saturates for C ≳ 75˚. The small η (≲ 0.1) and the tendency to form ribbons, which make the chains of the rods shorter in length, are the reasons for the absence of the extended clusters. The connectivity properties can also be studied by the analysis of the cluster size distribution n(s), where s is the size of the cluster, i.e. the number of rods belonging to the cluster. In Fig 7 we present the average cluster size distribution for different C with η = 0.2 [Fig 7(a)] and η = 0.4 [Fig 7(b)]. In general, n(s) decreases with increasing cluster size. The percolated configuration exhibits a single peak for large s, comparable to the system size, due to the finite size of the system considered in the simulations, and these states are denoted as random percolated [48,49]. For η = 0.2, the system is percolated for 45˚< C 60˚(see Fig 6(b)). Close to percolation (C ' 45˚) the n(s) curve presents a power law dependence, n(s) / s τ , with exponent τ ' −2.05, which is related to the well-known 2D random percolation prediction [50](τ = −187/91 ' −2.05) valid in the thermodynamic limit. Similar random percolated structures were found in Langevin Dynamics of functionalized colloids [51]. For η = 0.4, a similar s-dependence for n(s), also with τ ' −2.05, is found for C ' 30˚, where the rods aggregate together through head-to-tail bonds, similar as shown in Fig 1(c). In both, η = 0.2 and η = 0.4, a similar power law dependence n(s)*s τ is observed, but with different mechanism. For η = 0.2, the rods aggregate to each other through head-to-tail bonds and the chains of rods start to merge into one another in a random way. For η = 0.4, the rods are connected to each other as ribbon-like arrangements forming long chains, which, in turn, are randomly bonded to each other. The n(s) curves also confirm that the percolation transition for the system with larger C takes place for smaller η.
Another way to characterize the percolated structures is through the pair connectedness function g conn (r), defined as the conditional probability of finding a pair of particles separated by a distance r, both connected via a sequence of bonds, i.e., within the same cluster [52]. When an infinite cluster is present g conn (r) remains finite on every length scale. On the other hand, when a non-percolated structure is formed, we have g conn (r) ! 0 for finite distances. A theory of the pair connectedness function has been previously developed for fluids as well as for lattice systems when the presence of physical clusters of particles in the system is explicitly taken into account [53]. In the limit when the whole system forms a single cluster, the pair connectedness function matches the pair correlation function [Eq (11)].
In Fig 8(a) and 8(b) we show the pair connectedness function g conn (r) for packing fraction η = 0.2 and η = 0.4, respectively, and different values of C. An infinite-size percolated cluster is observed when g conn assumes non-zero values in the limit r ! 1, otherwise the cluster is not percolated. As can be observed the results are in agreement with those obtained through other distinct quantities, n(s) and S max .
In Fig 9 we briefly illustrate the effect of the increase of the aspect ratio on the percolation transition. For C = 15˚, larger values of l are associated with larger values of S max at η = 0.4. Specifically, the system for l = 3 at η = 0.4 is clearly percolated, whereas for l = 5 it is in the transition region, and for l = 7 the system is not percolated. This shows that the percolation transition is suppressed for larger values of aspect ratio. Such a result agrees with our previous study [19]. For C > 15˚, the behavior is the opposite, i.e., the percolation transition is enhanced for larger values of l. This suggests that when the head-to-tail arrangements are dominant, the attraction between the rods is weaker, as already discussed, and the geometrical effects resultant from the increase of the rods aspect ratio (in a 2D system) make it harder to form large clusters [19]. This changes when the ribbon-like arrangement comes into play, because this arrangement has a stronger attraction and the geometrical effects resultant from the larger η and larger l no longer hamper the formation of large clusters, since the side-by-side arrangements is expected to be more present in the higher packing fraction regime. Self-assembly of magnetic rods with tunable directional interaction

Effect of temperature
In this section we study the effect of temperature on the connectivity properties. The increase of temperature has a similar effect as to decrease the effective interparticle interaction. Although, from the perspective of reducing finite-size effects, it is advantageous to consider the percolation transition as a function of the density rather than as a function of temperature. However, it is worth to study the effect of temperature on the connectivity properties and to examine the stability of the clusters. It is expected that the percolation transition is hampered when temperature of an interacting system increases [52,54], since stronger thermal fluctuations may break the transient bonds.
In general, the ribbon-like arrangement is more stable against thermal fluctuations when compared to the head-to-tail structure. A specific example of the effect of temperature on the self-assembled configurations and on the percolation transition is presented in Fig 10(a), where the average size of the largest cluster (S max ) is shown for η = 0.4 and for three values of C. As pointed out previously, the system is percolated (S max ! 0.5) for all the three considered values of C for sufficiently small temperature (k B T/ = 0.1). In general, S max decreases with increasing k B T/, but how fast this happens depends on C. The lower the value of C, the lower the value of S max for a given temperature. As a consequence, we may conclude that the percolated configurations of rods with larger C are more stable against thermal fluctuations as compared to those with lower C, and this is due to the formation of ribbon-like arrangements (see Fig 10). When compared to the head-to-tail bond between two rods, the ribbon-like one is stronger, since all dipoles of one rod are more strongly attracted to the dipoles of the other rod for C > 45˚, forming local head-to-tail bonds. We present in Fig 10(b)-10(g) snapshots of the ribbon-like structures, with (C = 90˚and C = 60˚). For temperatures k B T/ = 0.1, 0.3, and 0.5, highlighting how the clusters melt when submitted to higher temperatures. From these results, we may conclude that the percolation transition is enhanced for a system of peapod-like rods with larger values of C even when it is analyzed as a function of temperature.

Conclusion
In summary, we investigated a two-dimensional system consisting of magnetic rods using Langevin Dynamics simulations. Each rod is composed of 3 soft beads having a central point-like dipole which interact via a DSS potential. This model is motivated by recent experimental [17] and theoretical [18,19] studies. A novelty of our system is the misalignment of the dipole moment of the individual beads with respect to the axial axis, opening the possibility of forming complex structures using nonspherical particles [12,[38][39][40]. Structural properties were investigated with particular attention to the dependence on the dipole's direction C and the packing fraction. We considered C ranging from C = 15˚to C = 90˚.
Due to the dipole-dipole head-to-tail assembly tendency, the increase of C produces a range of different configurations. As a consequence, we found that for C ! 45˚we observe ribbon-like structures which are energetically the most favorable. Given the preference of ribbonlike chain configurations we paid special attention to the appearance of a cluster extending over the whole system for a sufficient large packing fraction, i.e., the percolation configuration. We observed that the larger C, the stronger the attraction between rods, facilitating the clustering process and favoring the formation of the infinite cluster. As a consequence, the larger C, the smaller the packing fraction for the occurrence of the percolation transition. We found that by increasing the aspect ratio of the rods the packing fraction increases, and different geometrical effects are observed depending on the dominant kind of rod arrangement. For headto-tail arrangements, the percolation transition is hampered with increasing l. The opposite behavior was observed when the ribbon-like arrangement becomes dominant.
We also investigated the temperature dependence of the percolated configuration by analysing how the average size of the cluster depends on k B T/. We observed that percolated configurations with large C are thermodynamically more stable as a consequence of the formation of ribbon-like bonds, which is characterized by local head-to-tail arrangements.