Nanoparticle Shape Influences the Magnetic Response of Ferro-Colloids.

The interesting magnetic response of conventional ferro-colloid has proved extremely useful in a wide range of technical applications. Furthermore, the use of nano/micro- sized magnetic particles has proliferated cutting-edge medical research, such as drug targeting and hyperthermia. In order to diversify and improve the application of such systems, new avenues of functionality must be explored. Current efforts focus on incorporating directional interactions that are surplus to the intrinsic magnetic one. This additional directionality can be conveniently introduced by considering systems composed of magnetic particles of different shapes. Here we present a combined analytical and simulation study of permanently magnetized dipolar superball particles; a geometry that closely resembles magnetic cubes synthesized in experiments. We have focused on determining the initial magnetic susceptibility of these particles in dilute suspensions, seeking to quantify the effect of the superball shape parameter on the system response. In turn, we linked the computed susceptibilities to the system microstructure by analyzing cluster composition using a connectivity network analysis. Our study has shown that by increasing the shape parameter of these superball particles, one can alter the outcome of self-assembly processes, leading to the observation of an unanticipated decrease in the initial static magnetic susceptibility.

It is often the case that the significance of a newly realized material directly impinges on the range of capabilities it has. This statement simply comes down to a question of whether or not a certain degree of versatility is present. The ability to perform multiple tasks is crucial; simultaneous execution is certainly the pinnacle of such a goal, however it is often sufficient to simply switch between states that perform or behave differently. A quintessential example thereof is the ferro-colloid.[1] This magnetic liquid—a stable colloidal suspension of magnetic nanoparticles in a carrier liquid—represents a switchable system, regulated with the application of a magnetic field.[2] In addition, the capability to tune the exact state of the ferro-colloid through precise control of the applied magnetic field strength has only increased its utility.[3,4] The constituent magnetic particles, most commonly composed of some iron oxide variant, are undoubtedly the critical component characterizing the response of the fluid.

With this in mind, it is useful to make the link to a topic of significant current interest within colloidal science as a whole; specifically, the notion of short-range directional interactions that manifest in systems comprised of “shape-anisotropic” particles.[5,6] Of particular interest for us here is the distinctive phenomena arising in systems of magnetic particles that are themselves nonspherical.[7,8] The idea being that the nonspherically symmetric short-range interaction will be intrinsically coupled to the anisotropic magnetic interaction; in a sense, this scenario can be viewed as a competition between two sources of directionality. This competition will directly impact on the self-assembly mechanisms at play within the material.[9]

Given this fact, it is no surprise that a significant body of work on the realization of experimental magnetic particle systems of this type has occurred. These systems can generally be grouped into two distinct classes. Classification one: spherical symmetry of the particle geometry is typically retained, where some modification or symmetry breaking of the magnetic interaction is introduced; archetypical examples include Janus, capped, and patched colloids.[10−15] Classification two: particle geometry is the principal modification, and changes in the magnetic character are solely a consequence of this alteration; a wide variety of particle shapes have materialized, including cubes, rods, peanuts, and ellipsoids, to name but a few.[16−21] However, it is often the case that features of both appear in a single particle type, i.e., Janus magnetic rods and colloidal dumbbells with magnetic belts.[22,23]

The work presented here is concerned with the second of these three particle varieties, with our attention focused exclusively on cubic particles. Although we are only considering one particle type, an interesting additional consideration arises. Namely, the synthesis of colloidal magnetic cubes typically results in particles that are anything but perfectly cubic.[24] The curvature of the particles’ edges and vertices has a direct impact on the resulting behavior.[25] It follows that one has to explicitly account for this phenomenon when looking to accurately model these particles. Thankfully, a geometrical definition exists that has been shown to account for the particle shapes observed. This classification is termed a superball; the surface of which can be defined in a particle centered frame of reference as followsThe exact structure of a superball is determined by the value of the shape parameter q. Superballs with shape parameters lying in the domain q ∈ [1,∞] interpolate smoothly from a sphere (q = 1) to a perfect cube (q → ∞). The quantity h defines, independently of q, the size of the superball’s principal axes (particle height).

In recent years, a great deal about the behavior of superball particles has been elucidated. Numerical investigations that explored purely the effect of particle shape (nonmagnetic hard particles) highlighted interesting phase behavior.[26,27] Moreover, it should come as no surprise that predicting the influence of the shape parameter on the packing of these particles was among the first investigations to be conducted.[28] Experimental studies have not only addressed the crystal formation of superballs but also the rheological properties.[29−31] Predictions relating to the electrostatic and hydrodynamic characteristics of superball suspensions have also been reported.[32] The crystallization of magnetic superballs has also been investigated, with the caveat that the crystals were formed by sedimentation.[33] Furthermore, recent numerical studies on fluids of cube-like magnetic particles under the influence of an external field have begun to shed light on the implications of the particle shape on the aggregation present.[34,35]

One important aspect that has not been addressed to date in any detail is the explicit variation in magnetic response one would expect from spheres to cubes. Having noticed this omission, we chose to consider the following assertion: If we assume a superball ferro-colloid could be stabilized, then what would the differences in magnetic response be between it and a traditional ferro-colloid composed of spherical particles. The response of a magnetic material can be characterized by measuring the magnetization curve M(H): the variation of the magnitude of the magnetization M = |M|, as a function of the magnetic field strength H = |H|. In particular, the initial static magnetic susceptibility (ISMS) is a characteristic property of ferro-colloid; it is a measure of the magnetization curve within the linear response regime (i.e., small fields), χ0 = dM/dH|. Furthermore, the ISMS of magnetic fluids is a concept that is well studied not only from the point of view of experiment and theory but also computer simulation.[36] As such, it was useful to analyze this quantity for magnetic superball fluids to facilitate a meaningful comparison.

There are a number of theories that have been devised over the years that attempt to connect the macroscopic magnetic response of ferro-colloids and the interparticle interactions of the constituent particles. Langevin’s model of ideal (non-interacting) particles was the genesis of these efforts, where the ISMS for this ideal paramagnetic gas reads in terms of the number density of particles (ρ) as χ = 4πρλ/3.[37] This expression introduces an important quantity, the magnetic coupling parameter λ:that is, the energy per particle of two collinear dipoles (m = |m|), at close contact (h) relative to thermal energy, kT. This is the only parameter needed to systematically vary the strength of the magnetic interaction. The primary concern over the next century was how best to account for the interactions between the magnetic nanoparticles that Langevin had neglected. This question becomes even more pertinent with increasing λ. A selection of the most relevant attempts are summarized concisely by Kantorovich et al. and the references therein.[37]

Arguably the most successful of these theories to date is the so-called modified mean-field theory (MMFT); detailed comparisons between experimental data and simulation results, for not only monodisperse but also polydisperse systems, have confirmed its validity over the widest range of temperature and concentration.[37−39] The success of MMFT can perhaps be best described through its physical interpretation. MMFT provides a framework to account for the effective field acting on a particular particle; the effective field depends not only on the applied external field but also on all of the other magnetic particles in the system. The expression for the ISMS found following a rigorous derivation from first-principles reads as follows:it being a second-order correction to the model proposed by Langevin.[39] The initial statement of MMFT for the ISMS up to first order actually came about as an empirical modification to the Weiss model.[37,38] Subsequent developments of MMFT incorporated the important effects of chain formation into this framework.[40] It follows that MMFT, including chain formation, is a prime candidate to apply to fluids of magnetic superballs, albeit with a few adaptive modifications.

In order to address the lack of insight in this area, we present here numerical investigations exploring the magnetic response of relatively low-density fluids of magnetic superballs. The field produced by a magnetic superball was modeled using a single dipole placed at the particle’s center of mass and oriented in the [001] crystallographic direction. The study of colloidal magnetic suspensions often focuses on tailoring the system’s magnetic response. This procedure can involve targeting particular behavior under certain physical conditions or, as the case here, attempting to realize systems that react stronger to an external magnetic field. Cube-like particles were considered possible contenders for the realization of such an aim, a conjecture that was based on results from previous work which reported that the preference for magnetically “reactive” chain aggregates in the ground state of superball particles increases with growing q.[25,41] In other words, the nonzero dipole moment of the observed chain structures could, in principle, lead to a more field-susceptible system. Furthermore, there is a real prospect that the documented closing transition for spherical particles into ring aggregates, which induces a abrupt decrease in ISMS with reducing temperature, has the possibility of not occurring for cubic particles.[42] It transpires that the thermal stability of the structural motifs known to constitute the ground state of superball particles is a key consideration at finite temperature, an effect which we quantify in this work.[25] The crux of our investigation became the following: How do modest increases in the value of the superball shape parameter (q = 1.0, 1.25, 1.50) impact the resulting ISMS of the corresponding magnetic suspensions?

We have also sought to adapt MMFT, with the aim of accounting for the specific shape parameter of the superballs in question. Moreover, we have studied the microstructure of each suspension via an analysis of the cluster composition, allowing us to reconcile the effect of cluster formation on the calculated macroscopic magnetic response. Even though the use of the magnetic dipole–dipole potential to characterize the long-range interactions of nanocubes seems rather bold, it can still be considered valid for either single domain magnetic spheres embedded within nonmagnetic cubic shells[20] or can be regarded as the leading term for explicitly cubic nanocolloids with remanent magnetization.[43]

It is important to acknowledge at this point that particle polydispersity can have a sizable effect on the magnetic response of dipolar systems, a fact that is exacerbated in this case as there are two sources of polydispersity: size and shape. Due to the significant complications arising in the multicomponent system that would be required to model these effects, we have limited our study to a strictly monodisperse regime. Nevertheless, it is known from studies of polydisperse dipolar hard spheres that a poisoning effect occurs, which inhibits cluster formation and thus magnetic response; the same scenario would likely manifest itself in this case.[44] The effect of shape polydispersity, on the other hand, is somewhat harder to predict. What is clear, however, is that the magnetic response of these systems will strongly depend on the exact distribution function characterizing the variation in q.

The results of our studies are found in the subsequent section. We present the simulation results relating to the ISMS of the suspensions, including a direct comparison to the predictions of an adapted MMFT. Appearing alongside is a description of the cluster analysis algorithm and the resulting cluster distributions. We conclude our discussion of this work with a summary of our findings, highlighting the particular significance the results could have for an experimentally realized version of our system. The model and methods we have used are surmised at the end of our contribution, where the important features of the established particle model and simulation technique used are detailed, accompanied by comments on the alterations required to adapt these to calculations of bulk properties.

Results and Discussion

Initial Static Magnetic Susceptibility

The correlations between dipolar particles are particularly important for theoretical considerations relating to any kind of measure of magnetic response. Moreover, the observation of chain formation in both experiment and simulation emphasizes the need to account for aggregation in the system.[36,45] Increases in the ISMS observed are a direct result of the highly correlated particles found in chain aggregates. Mendelev and Ivanov extended the formulation of MMFT to account in part for the formation of chains.[40] The expression for the resulting ISMS was formulated asby introducing parameters describing both a measure of the probability of bonding p and two-particle correlations K. The first-order variant of eq is the prefactor of the chain forming ratio. The derivation of MMFT assumes the magnetic particles to be spherical.

To account for the explicit shape of superball particles, we have devised a simple manipulation to this expression. Namely, we proposed that the bonding probability, p(q), must be dependent on the value of the shape parameter. Within the established scheme, the bonding probability arises as the Lagrange multiplier required for the minimization of the system free energy density in the presence of a mass balance condition.[37] Its functional form depends on the system volume fraction φ(q) and the two-particle (dimer) partition function , which consequently are now both shape parameter dependent:The system volume fraction is expressed in familiar terms for superballs as follows:where N is the number particles and V the system volume, νsb(q) is the dimensionless single particle superball volume, and Γ(x) is the gamma function of argument x. For the dimer partition function we use the following expression:The dimer partition function for spheres is equivalent to that given in ref (40). The numerical factor f(q) appearing is the result of a linear interpolation from spheres to cubes in terms of particle volume. In contrast, we retain the q-independent definition of the correlation coefficient from ref (40):For those interested in the details of our modification please consult the Supporting Information.

Within the theoretical framework, we have just outlined that the key particle-shape-dependent quantity is of course the dimer partition function. As such, it is useful to visualize its behavior and note the differences resulting from changes in the value of q. To this end, the dimer partition functions for the q values under consideration are found in the upper plot of Figure . In addition, the corresponding bonding probabilities appear underneath and are plotted for systems with a reduced number density of ρ* = 0.075, illustrating the effect of q on a more physical variable. Due to the modest perturbations from sphericity considered here, we have chosen to show one further shape parameter (q = 5.00) simply for illustrative purposes. One can see that the overall trend is for to increase in line with increasing λ. Importantly, however, one also notices that as q is increased, the number of accessible states drops ( decreases) for all values of the magnetic coupling parameter. Turning to the bonding probability, we observe the same pattern; p increases with increasing λ, but falls when q is increased for fixed λ. The expected convergence of p to unity is observed with increasing λ and for each value of q. The function and in turn p directly determine the influence of chain aggregation in the calculation of the ISMS. Already, the variation of p and as a function of shape parameter predicts and in fact induces a fall in ISMS.

Figure 1

Plots of the analytically derived expressions for the dimer partition function (upper plot) and bonding probability p(q) (lower plot). The curves are shown for the shape parameters under investigation, q = 1.00, 1.25, 1.50, in addition to the illustrative supplementary value of q = 5.00. The coloring of the curves follows the color bar appearing at the top of the figure. Please note that a log scale has been used on the ordinate axis.

From the perspective of our simulations, we compute numerically the ISMS on the basis of the fluctuation–dissipation theorem. We have used the following fluctuation formula:The reduced number density of particles in the system is given by ρ* = ρh3, and the instantaneous magnetization of all dipole moments is determined accordingly, M = ∑m. The upmost care must be taken in the computation of the magnetic fluctuations required in eq . The dynamics of dipolar systems, especially for strongly interacting particles (high λ and/or high ρ*), can be extremely slow. Thus, simulations typically ran for significant periods, enabling the relaxation of the system and the equilibration of these fluctuations. The length of this “significant period” was determined on a case by case basis, whereby we ensured that the residual magnetization |M| in the system becomes negligible within the limits of statistical accuracy.

The parameter space of our investigation was based on the following considerations. First, we were most interested in what effect small perturbations from spherical particle geometry would have. That is why we only considered relatively small increases in the shape parameter, specifically q = (1.25, 1.50). Furthermore, in the interest of an effective comparison, we performed the same set of simulations for spheres (q = 1.00). Second, as previously stated, we have considered these magnetic superball suspensions in a low-density regime with ρ* = (0.01–0.1) in increments of 0.025; in addition, a pair of higher values ρ* = (0.15, 0.2) were also considered, to ascertain the effect of higher concentrations. Finally, magnetic coupling parameters, ranging λ = 1–5 in increments of 1, were considered in order to navigate a large spectrum of self-assembly scenarios: from negligible via transient to significant system aggregation. The full details of the model and methods we have used in our simulations can be found in the Model and Methods section.

To begin, it is beneficial to look at the variation in ISMS as a function of our three basic state variables (q, ρ*, λ). However, according to the superball definition, the particle height h remains constant when changing q, yet the volume of a single particle does not. Therefore, after some consideration it makes sense to map the data from number density to volume fraction, in order to affect a more meaningful comparison. The mapping is a simple change of variable according to eq . This conversion allows us to collate all the data into a single detailed plot of χ0vs φ that appears in the main portion of Figure . The points represent data from simulation, and the curves are the theoretical predictions of our slightly modified version of the chain formation variant from eq . The legend at the top of Figure shows the allocations of the line and symbol types for each q. The color coding of both the data points and curves relates to the value of the magnetic coupling parameter, according to the color bar shown. A logarithmic scale has been used on the ordinate axis to facilitate a clearer distinction between the data points and curves shown. The error bars appearing are a measure of the standard deviation of the mean (please see the Model and Methods section).

Figure 2

Plot of χ0vs φ. Symbols indicate data attained from simulations, and lines denote the theoretical predictions. The line and symbol identifications obey the legend appearing at the top of the figure. The color code relates to the value of the magnetic coupling parameter, according to the color bar at the top of the figure. The error bars shown denote the standard deviation of the mean. Inset: plot of Δχ0vs ρ*. Symbols and coloring are as described for the main plot, however the lines appearing in the inset are purely to aid the eye. The error bars shown have been propagated from the corresponding errors in the main plot.

The general trend observed for each q was of increasing χ0 with growing λ for fixed φ, as to be expected. The most striking feature of this plot is the observed decrease in χ0, for fixed λ, as q is increased; an observation evident for all values of λ. In terms of the theoretical expression, it performs exceptionally well in the low coupling regime (λ ∼ 1, 2) even to the highest φ considered. One also notices as we move to the higher coupling regime (λ ∼ 3, 4), the agreement is still very good for each value of q, albeit for the deviations that begin to occur as φ reaches its highest values; interestingly, however, the accuracy improves as the superballs become more cubic. For the highest coupling λ = 5, a significant difference between the simulated and calculated data develops (except for the lowest φ), which is in line with the assumptions/limitations of the chain formation model in general.

To further access the differences in static response that result from increasingly cubic particle structure, we have considered the following measure:the difference between the ISMS of spheres and cubic particles, q = (1.25, 1.50). This measure is only illustrative at the higher values of λ considered, once the geometry of the particle begins to properly influence the self-assembly that occurs. With this in mind, we have plotted Δχ0(q) for λ = (3–5) in the inset of Figure ; the functions are plotted in terms of reduced density to allow the differences to be calculated in the first place. The plot exemplifies the changes in system response that arise due to the variations in particle geometry; changes, which are exacerbated at higher values of λ, especially in the lower half of the density range. Another point of interest that this measure raises is the possibility of a regime change in the evolution of the ISMS as ρ* or φ increases. Namely, one can see for λ = (4, 5) in particular, we have a clear maximum in Δχ0. This suggests that after a certain value of ρ* (depending on λ), the ISMS of cube-like particles seems to be recovering to values more comparable to that of spherical particles. This posited regime change is also evident in the main plot of Figure , where for λ = (3, 4) we can see the value of φ where the chain formation model begins to abruptly fail. It seems then that the chain formation model can accurately describe the particle correlations that manifest in the initial mechanism, but fails to do so once this secondary state has been entered. We believe this second region, the onset of which is dependent on all three system parameters (q, ρ*, λ), results from the evolving nature of the agglomeration occurring as the systems become increasingly concentrated and correlated; perhaps ring formation is beginning to exert some influence.

One further mapping, to the Langevin susceptibility χ, allows for another commensurate comparison. In particular, it emphasizes variations in particle correlations. This plot of χ0vs χ appears in the upper portion of Figure ; the line, symbol and color identifications are exactly equivalent to those described for Figure . The dashed-dotted black line gives the linear Langevin ISMS law, which one should recall neglects the effect of magnetic correlations between particles. We have limited the range to 0 < χ < 2.0 so that the densely packed region at low χ is easier to see. With this in mind, one can see that for values of χ < 0.5, the linear relationship holds, indicating the absence of strongly correlated particles. For χ > 1.0, we begin to see large discrepancies between the data and the linear curve, indicating the increasing importance of the magnetic correlations. Furthermore, one can state that for a given χ (at higher values of λ), the same trend of lowering χ0 for increasing q is observed, unmistakably showing that systems at the equivalent state-point are less correlated for superballs that are more cubic. One should note the fact that the theoretical curves for each q at the various λ show little difference. This is due to the fact that there is no explicit variation in how the correlations are calculated as q is increased.

Figure 3

Upper plot: χ0vs χ. The color bar above the plot indicates the value of magnetic coupling parameter. The black dashed-dotted line represents the linear Langevin ISMS law. Lower plot: χ0vs λ for equal νsb, where λ has been scaled by a factor of 1/(h*)3, which depends on the value of q. The color bar at the top left of the plot denotes the respective system density. Please note the use of a log scale on the ordinate axis of this plot. The line and symbol identifications obey the legends appearing in each plot, respectively. The error bars in both plots denote the standard deviation of the mean.

The approach we have used here was chosen as such to scale in terms of comparable interaction energies, rather than comparable particle volumes. Yet in experiment, particle volume and the magnetic moment are directly related, m ∝ ν. Thus, in general, λ can not be kept fixed as particle volume increases. To address this issue, we performed a number of complementary simulations for systems at ρ* = (0.05, 0.075), for fixed values of the magnetic moment (m*)2 = (1–5) in integer steps. We fixed the particle volume for each value of q to νsb = π/6, which consequently varies the length of the principal axis, h* = (1.0, 0.9382, 0.9026) for q = (1.0,1.25,150), respectively. The variation of the ISMS with λ for these systems is shown in the lower plot of Figure , where the symbol and line identifications are given by the legend in the lower right-hand corner, and each density is colored according to the inset color bar appearing at the top left. The ordinate axis of this plot uses a log scale. It is crucial to note that as the size of h* now varies as a function of q, λ should adjust accordingly, i.e., by a factor of 1/(h*)3. We can see from the figure in question further evidence of the decrease in susceptibility, using the appropriately scaled values of λ. The theoretical predictions again do well, but also begin to break down as the secondary regime is entered at λ ∼ 4.

Interestingly, however, this does not tell the full story. One further observation, at constant density, can be made. It is evident from the figure that systems with low m* have an ISMS that does not significantly change with increasing q, while for higher values of m*, the ISMS is seen to increase as a function of q. This implies the following extremely interesting assertion. Suppose one could make a number of stable colloidal suspensions in the laboratory using superball particles of precisely equal volumes with different shape parameters. If one were to make no assumptions regarding the strength of the interactions involved, the measured ISMS of each suspension would actually be the same, if not higher, as the particles become more cubic. At first glance, this seems to be a major contradiction. However, this comparison is not made on an equal basis in terms of the physics of the system we have modeled; namely, we consider a colloidal particle with a point dipole placed in its center. Therefore, the increase measured, rather than originating from the effect of the particle geometry, would simply be related to the fact that, at equal particle volumes, the length of the particle’s principle axis decreases as q increases. Hence, the observation is simply due to the fact that the dipoles of each particle can get closer to each other as they become more cubic. Thus, to properly assess the effect of the particle geometry, a comparison in terms of interaction energy is the correct approach. As such, the rest of the manuscript focuses on the analysis of systems based on this foundation, i.e., particles whose principal axes are of equal lengths, h* = 1.

It is useful at this point to provide a brief interlude to summarize what we have elucidated so far. Primarily, we now know that the ISMS of a bulk suspension lowers as the constituent magnetic superballs become more cube-like; a fact that is neatly encapsulated by our basic modification of the chain formation model for spheres; a framework, which we show is valid for a preliminary regime encompassing, in general, the lower end of the range of φ and λ ⩽ 4. Second, the onset of what we believe is a subsequent regime has been identified, but as yet not fully explained. Third, the observation of evidently higher correlations for more spherical superballs suggests that the same trend should be present in the structure of the system aggregation one should observe. To clarify a number of these statements, we chose to ascertain a more detailed picture of the microstructure present, allowing us to rationalize and account for the observed macroscopic magnetic response.

Cluster Analysis and Microstructure

In order to provide some kind of explanation for the drop in ISMS we have observed, we constructed connectivity networks representative of the system evolution at each state point. Thus, a detailed understanding of the cluster formation could be achieved. The first application of this connectivity network method to dipolar systems was in relation to polymeric style magnetic brushes composed of dipolar filaments.[46] The methodology applied here follows closely this implementation, and those looking for the intricacies and full terminology of the procedure should consult ref (46) in detail. Furthermore, we recently trialled the applicability of this approach in relation to suspensions of magnetic spheres and cubes.[47] This work provided the proof of concept necessary to allow the successful application of this technique to the current concern.

The algorithm used here follows closely that described in ref (47), and what follows is a number of reminders regarding definitions and bonding criteria purely for the reader’s convenience. Within a given system configuration, each particle forms a single vertex of a connectivity graph. If two particles are bonded to one another, an edge is drawn between the two corresponding vertices. Two particles are said to be bonded if the following criteria are met: the two particles must be within a certain distance of one another, i.e. ⩽1.1, and the dipolar interaction energy between them must be favorable, i.e. < 0. The distance criterion is derived from the two particle interaction potential; it represents the distance at which two collinear dipoles are at 75% of the attractive well-depth. Once all the edges have been mapped out, the connected subgraphs constitute the clusters that have formed in the system, indicating immediately their size and general topology. A particularly useful measure of interest is termed the vertex degree, it being the number of edges (bonds) associated with a particular vertex (particle). As such, the vertex degree can give us an excellent idea about the propensity with which particles like to form bonds. As with any system observable, all the measurements taken regarding cluster composition are reported as ensemble averages. It should be noted, as periodic boundary conditions were in effect during the simulations, extra care was required when identifying bonded particles. Specifically, clusters were allowed to form across the periodic boundaries of the simulation box, ensuring that no clusters were artificially truncated.

As an initial measure of the aggregation occurring in these magnetic superball systems, we will turn to a simple pictorial assessment, namely, the series of simulation snapshots appearing in Figure ; the snapshot array is for the set of systems with ρ* = 0.075. Each row of the array is for a single q, which increases in a downward direction, and each column is for a single λ increasing from left to right. Each snapshot shown is the final configuration recorded from the production period of the respective simulation runs. The coloring of the superballs is a result of the cluster analysis we have performed. A particle is colored according to the cluster to which it belongs; the cluster color is determined by the color bar appearing on the right-hand side of the figure, noting that monomers are colored white. This way we can get an immediate impression of the amount and scale of the agglomeration at each state-point. Along a single row, we see the aggregation increase for each value of q in line with the growing value of λ in terms of not only the number of clusters but also the size of individual clusters. Down a column, however, the opposite is true; the abundance and length of clusters reduces as q increases for fixed λ. Considering the snapshot array in its entirety, the propensity for cluster formation is in general much higher for superballs that are more spherically shaped. In the interest of completeness, the equivalent snapshot arrays for the other values of ρ* can be found in the Supporting Information. The observations made bode well for the links that will be drawn to the magnetic response later on; beforehand, more quantitative measures of the bonding are required.

Figure 4

A snapshot array of state points with a density of ρ* = 0.075. Each row corresponds to state points with the annotated value of q, and each column corresponds to the annotated value of λ. The coloring of clusters is in accordance with the color bar appearing on the right-hand side of the figure, where monomers are colored white.

Two further quantitative measures that illustrate the variation in cluster formation across the parameter space: the average cluster size ⟨N⟩ (excluding monomers) and the aggregation percentage (Agg–%), i.e., the percentage of particles in the system forming at least one bond. The first measure is of course an indicator of the size that clusters in the system reach, while the second is an indicator of the overall connectivity within the suspension. These two quantities are displayed in the form of bar charts appearing in Figure . Both bar charts are organized in the same way: the quantity of interest appears on the ordinate axis; the abscissa organizes the bars into seven groups corresponding to the value of ρ* indicated; within each group the value of q is indicated by the bar hatching, circles (q = 1.00), horizontal lines (q = 1.25), and diagonal lines (q = 1.50); the value of λ is given according to the color bar appearing on the right of the figure; finally, the extra spacing appearing in the groupings for ρ* = 0.15, 0.2 is due to the absence of λ = 5 for these densities. Presenting the data in this fashion allows one to identify the trends across each variable ρ*, q, and λ.

Figure 5

Top: bar chart of ⟨N⟩ vs ρ*. Bottom: bar chart of Agg–% vs ρ*. The respective density value of each bar chart set is given at the top of each plot. The hatching of the bars indicates the value of the shape parameter according to the following identifications: circles, q = 1.00; horizontal lines, q = 1.25; and diagonal lines, q = 1.50. The color coding for the respective λ in these plots obeys the color bar appearing at the right-hand side.

Beginning with average cluster size, we observe for a given ρ* and q the increase of ⟨N⟩ for increasing λ, as one would expect due to the increasing strength of the magnetic interaction. For a given ρ*, the tendency for ⟨N⟩ to decrease with increasing q is clearly seen for all values of λ. Furthermore, we also observe a slight (if not almost negligible) increase in ⟨N⟩ as a function of ρ* at low values of λ for each q, whereas the rate of increase accelerates at higher values of λ. Moving to the aggregation percentage we can see that exactly the same trends are prevalent, although much easier to elucidate, as the variation of this quantity between state points is much stronger. The trends we have identified corroborate our visual inspection of the snapshots from Figure (and in the Supporting Information), albeit the snapshots only give us an instantaneous interpretation of the system in comparison to these time-averaged measures. Regardless, the equivalent behavior is evident from both: the onset of agglomeration in magnetic superball suspensions is delayed to higher values of λ as the shape parameter of the particle is increased, and for a given ρ* and λ, the connectivity of the suspension and the size of clusters formed decreases as the particles become more cubic.

We now turn to the vertex degree δ that was alluded to earlier. It maps out clearly the bonding present in the systems. We have chosen to present the data using histograms in order to give an idea of the average distribution of δ at each state point. The histograms are shown in Figure . Each individual plot is for the denoted pair of values (q, ρ*), where ρ* increases downward in three plot intervals. Within a single plot is a histogram for each value of λ, once more colored according to the color bar at the top of the figure; an equivalently colored vertical line denotes the average value of the corresponding distribution. The ordinate axis represents the probability of a particular value of δ, where the abscissa is the value of δ binned according to the Freedman–Diaconis rule. The vertex degree was recorded and averaged for each particle in a given simulation. Therefore, a single histogram for given (λ, q, ρ*) represents the distribution of the particles’ average δ.

Figure 6

Plots of δ histograms. The parameter indicator appearing at the top right of each plot denotes the values (q, ρ*). The color bar at the top of the figure indicates the value of magnetic coupling parameter. The vertical colored lines denote the average value of the corresponding distribution.

Immediately, one can see from each plot the progression of the distribution to higher values of δ, as λ is increased. Furthermore, for each individual value of ρ* the extent of this progression is hindered with increasing q. Interestingly, the progression for each set of fixed λ and q appears to begin to saturate as a function of ρ*. This array of plots illustrates definitively the influence of q, specifically on the bonding taking place in each system. As we can interpret δ as the number of bonds a single superball likes to form, it is explicitly clear from the histograms that the increase of q significantly hampers bond formation. Additionally, this analysis allows us to make some general comments on the probable topologies of the clusters formed. In particular, the data suggest the tendency for the formation of nonbranched cluster types i.e., chains and rings. The values of δ measured are all significantly <2, suggesting the probability of a branched structure is extremely low, owing to the requirement that a superball particle acting as a branching point would need δ ⩾ 3. Furthermore, state points with δ ∼ 1.0 will actually tend to be abundant in dimers; the value of δ implies a reluctance of particles to form subsequent bonds.

To conclude our analysis we looked inside each bond at the configuration of the constituent dipoles. In particular, we analyzed the distribution of the interdipole angle θ and both interdipole–displacement angles θ(1) and θ(2). These angles are defined for two bonded particles 1 and 2 asEach angle runs over a domain of (0, π). No constraints were imposed upon these angles, avoiding any unintentional bias being applied to the direction of the displacement vector. Its direction was treated in the same manner throughout the analysis, i.e., always from particle 1 to particle 2 of the bond, irrespective of the dipole orientations. Calculation of these angles not only provides an indication of the symmetry of clusters but also visualizes how the bonds fluctuate. This wealth of information can be attained by constructing bivariate histograms for pairs of these angles. Specifically, we have plotted the probability density functions, f(1) = f(θ,θ(1)) and f = f(θ(1),θ(2)), binned according to a two-dimensional variant of the bin optimization routine outlined by Shimazaki and Shinomoto.[48] Due to symmetry, the probability density function f(2) = f(θ,θ(2)) is qualitatively equivalent to f(1) and thus provides no further insight. The histograms for f(1) and f are plotted in Figures and 8, respectively. These plot-arrays appear for the same parameter set used in Figure , i.e., ρ*=0.075, with q increasing down each column and λ increasing along each row. Each bin is colored according to the corresponding value of the probability density function; the color bars used are shown on the right of each figure. As the propensity for bonding increases with λ and decreases with q, the bin size adjusts accordingly. Quantities such as these are not strong functions of the system density, and as such, no further information or understanding is gained from the presentation of systems at the other densities considered.

Figure 7

Histogram array of the probability density function f(1), of the angle pair (θ,θ1). The state points displayed have a density of ρ* = 0.075. Each row corresponds to state points with the annotated value of q, and each column corresponds to the annotated value of λ. The coloring of individual bins is in accordance with the color bar appearing on the right-hand side of the figure and represent the value of f(1) at that particular angle combination.

Figure 8

Histogram array of the probability density function f, of the angle pair (θ(1),θ(2)). The state points displayed have a density of ρ* = 0.075. Each row corresponds to state points with the annotated value of q, and each column corresponds to the annotated value of λ. The coloring of individual bins is in accordance with the color bar appearing on the right-hand side of the figure and represents the value of f at that particular angle combination.

The evolution of these histograms as a function of λ and q tells us a great deal about the system. If we consider the first pair of angles (θ,θ1), shown in Figure , we see the emergence of two peaks as λ is increased, an observation that occurs for each value of q. The sharpness and height of these peaks strongly increase as a function of q. Furthermore, we see regions of the parameter space become void of bonds with increasing q, particularly for higher values of λ. Focusing on the three histograms for λ = 5, we can see how the peaks proceed to smaller values of θ, indicating that the dipoles in the bonds become better aligned with one another. The peaks also proceed to the extremities of the θ(2) range (i.e., 0 and π), which suggests that the dipoles within the bond become increasingly colinear. The reason for the appearance of two different peak positions in θ(2) is simply due to scenarios where the displacement vector is treated in the opposite direction to that of the dipoles, resulting in angles also in the vicinity of π, not just 0. This fact also provides an additional equilibration check, whereby the distribution is extremely symmetric across the line θ(2) = π/2, suggesting no remanent magnetization is present in the system. In general, these histogram arrays suggest that the diffuse nature of the peaks for low q should allow for not only straight linear cluster formation but also for the presence of flexible configurations. As q increases, the dominance of straighter chains is cemented, because the dipole configurations within bonds favor parallel colinear confirmations.

Moving to the second angle pair, in Figure , we can make some further interesting observations. The same generic trends observed for the peak formation in the first angle pair also hold here. We have two peaks developing with increasing λ, which sharpen and tend toward (0,0) and (π,π) as q increases. This confirms the increasingly straight shape of chains as q increases. The appearance of bond-deficient regions of the parameter space is again observed, but to a lesser extent than for the previous combination. This indicates that increasing q acts to restrict the mutual orientation of the dipoles (θ) to a far greater extent than the displacement related angles. This validates, in part, our assumptions regarding the calculation of the dimer partition function, where we imposed this very criterion. A perfectly straight chain can only be realized by constraining all three angles. This is crucial when considering, in addition to the head-to-tail dipole configuration, the antiparallel pair configuration that dipoles can adopt. One can access their prevalence in these systems by analyzing both sets of histograms in the proximity of (θ,θ1) = (π,π/2) and (θ1,θ2) = (π/2,π/2). In doing so, it becomes clear that this region of the parameter space is certainly not unpopulated. The number of bonds adopting configurations close to an antiparallel pair, however, is significantly lower that of bonds adopting head-to-tail configurations. This disparity grows as q increases. Therefore, we can safely say that the influence of antiparallel dipoles greatly diminishes as the dominance of straights chains establishes itself with increasing q.

At this stage, we hope that the parallels between the ISMS of each state point and the corresponding microstructure we elucidated there are clearly discernible. It is evident that each of the measures discussed are self-consistent and substantiate the observed trend in the ISMS. The prevalent notion to be gained from ⟨N⟩ and Agg–% in Figure is for a given ρ*, both quantities decrease with increasing q, across each value of λ considered. This highlights the fact that the self-assembly taking place is inhibited as the colloidal particles become more cubic. These two facts themselves provide the bulk of the explanation for the decrease in ISMS. As the size and concentration of clusters decrease, fewer fluctuations in the instantaneous magnetization occur, which, as we are in the linear response regime, must result in a lowering of the ISMS. The reduction in the average number of bonds formed for more cubic particles, as exemplified in Figure , provides further evidence for the reducing probability of bonds forming. Moreover, Figures and 8 illustrate the fact that fewer low-energy (less negative) dipole configurations are available as q increases. In one sense, this can be seen from the perspective of the cubic particles as a delay in the onset of self-assembly to ever higher values of λ. In other words, the more cubic a particle is, the stronger the magnetic interaction required for bonding to occur. This assertion explains, in part, why the theoretical framework performs objectively better the more cubic the particles become; it is simply the case that there are less correlations in the system, which means the ones present are easier for the framework to encapsulate.

One question remains however, namely, what is the reason for the apparent aversion to bonding that cubic particles have in comparison to spheres at the same λ? The appropriate reasoning can be found by considering the concept of bonding volume. This bonding volume relates to the proportion of phase space where the formation of a bond between two superballs can take place. In this context, this is the region where the separation and orientation of two superball dipoles is such that it is sufficiently favorable to hold the two particles together. It follows that this bonding volume changes as a function of the shape parameter; it decreases and then saturates as the particle tends toward a perfect cube. This reduction in the available phase space means that it becomes immediately more difficult for the particles to form bonds, even for these moderate perturbations from sphericity. The presence of this effect is evident from the histograms in Figures and 8. It is clear that the number of angle combinations that bonds can adopt reduces significantly as q increases. This is a physical demonstration of the reduction in accessible phase space volume. These data reinforce the idea of decreasing bonding probability with growing q, primarily due to the number of suitable stable configurations decreasing.

At the beginning of this manuscript we posited the prospect of cubic particles creating more field susceptible suspensions, by way of their increased preference for chain aggregates. Not only have we shown this preconception to be false, we have provided the reasons why this is the case. In particular, even although increasing q shifts the cluster equilibrium from flux-closure rings to field-sensitive chains, we have shown that the probability of chain formation actually decreases as q increases. This observation comes with the caveat that under the physical conditions considered, the explicit competition between chains and rings is expected to be low. However, the decreasing difference between the ISMS of cube-like and spherical particles at λ ∼ 4–5 suggests that, as the precursor to ring formation begins for spheres, the evidently straighter chain clusters of cubic particles prevent the mechanism from occurring to the same extent. This presents the opportunity for the ISMS of cubic particle systems to, in a sense, catch-up. In the working region of parameters considered here, the limitations in effect for cubic particles delay the onset of cluster formation, resulting in objectively less correlation and thus the lower susceptibility we observed.

Conclusion

The rationale behind this work was a simple one: What effect does particle shape, cubic geometry in particular, have on the corresponding magnetic response of a stable magnetic fluid of colloidal superballs? The desire to address this question was motivated by the suspicion that owing to the preference of cubic particles with a [001] dipole to form chain aggregates, higher values of the ISMS could have been achievable. A supposition based on the fact that chains are a much more field-favorable cluster topology than rings.[41] At the “ambient” temperature conditions considered here, the contribution of ring structures to the susceptibility of spheres is somewhat negligible. This shifts the dependency of the ISMS from the explicit topology of clusters that result from self-assembly to the two-particle bonding process itself. In this way, we managed to attribute the q-induced drop in ISMS to this basic consideration, correcting in the process our initial preconceptions.

Our understanding of the system was enriched by a detailed analysis of the microstructure at each state point. This analysis not only corroborated the conclusions drawn from the computation of the ISMS but also facilitated the understanding of the self-assembly mechanisms in effect. The ability to visualize the aggregation in the system allowed for the significantly less correlated cubic systems to be easily identified. Furthermore, our modifications to the existing MMFT have proved valid. The resulting framework—particularly at low volume fractions and moderate magnetic coupling—faithfully represents the variation in the response of the system as a function of the shape parameter.

At first glance, it would appear that the tempered response of more cube-like particles is somewhat of a set back toward the aim of producing more field-responsive magnetic liquids. Be that as it may, there are a number of scenarios we can foresee where the opportunity for cubic particles to exhibit an advantage over their spherical counterparts may persist. The indicators of such prospects are found in the results of this work, but further study is required to confirm these conjectures. The first of these scenarios has been alluded to previously, namely, the effect of temperature. Although the persistence of ring structures for similarly q-valued superballs was recently reported, the likelihood of the abrupt decrease in ISMS with decreasing temperature is unlikely.[25] Thus, the prospect of cubic particles being more suitable for low temperature conditions is still very much in the cards. The second scenario that comes to mind is one where a sharper transition in terms of magnetic coupling strength is required. Careful consultation of the cluster analysis (numerical and pictorial) seems to indicate the possibility that the transition from a dispersed to an aggregated system occurs more quickly (i.e., over a smaller range of λ) as the particles become more cubic. The evolution in the aggregation of the spherical system appears more gradual. Thus, one can envisage the situation where a fluid with a sharper transition would benefit applications where the switchable nature of a ferro-colloid is of primary concern. Two state binary operation may allow for greater precision to be achieved.

Model and Methods

Particle Model

In our simulations we have approximated the explicit shape of superball particles using a scheme developed in ref (25); an extension to a preliminary model described in ref (41). The model in question utilizes constituent spheres to construct the superball surface according to eq . This is successfully achieved by placing the spheres at significant locations within the superball geometry, namely the vertices and the midpoint of each edge. The radii of the spheres are then set to match the radius of curvature of the superball at the given location. In this manner, one can mimic the superball geometry for low values of the shape parameter. As we have only considered modest increases in q, this model is more than suitable for our purposes. We have used a total of 21 spheres to construct our superball particles (8 vertices and 12 edges; plus one central sphere of diameter h*).

There are two interaction types present in our model: steric and magnetic. The nonmagnetic interaction between two superballs is simply that of excluded volume. For the purposes of our analytical calculations, the particles are treated as hard bodies. This scenario is emulated in simulation using Weeks–Chandler–Anderson potentials that act between the spherical sites of opposing superballs.[49] As discussed in ref (25), the cutoff radius r for the steric potentials is defined as r = 21/6σ; the interaction parameter ϵ, specific to the short-range potentials, was set to 1000; and all of the average site diameters σ were scaled by 2–1/6. This increases the similarity to a perfect hard particle scenario by steepening the potential and shifting the effective hard particle cutoff to unity. Therefore, the particles can be effectively considered as hard particles of size h*. This type of short-range interaction is a straightforward representation of the steric or ionic stabilization required in experiment to produce a colloidal suspension. A commonly asked question relates to our use of this model as apposed to an analytical overlap approach. The answer is discussed in some detail in ref (25); the arguments relate essentially to the efficiency of calculating the magnetic interaction, whereby advanced algorithms and parallelization can be employed in a more advantageous manner. Furthermore, it is important to recall that the hard/softness of steric potentials do not impact greatly on the resulting magnetic character of dipolar systems, especially for the range of magnetic coupling parameter we have considered.[39]

The magnetic character of the particles is described using a dipolar approximation. A single dipole was placed at the center of mass of each superball and oriented to lie along the [001] crystallographic axis, i.e., from face-to-face. Thus, the magnetic interaction between two dipolar superballs was governed by the familiar dipolar formulation. This choice of magnetic configuration makes a number of key assumptions about the particles in relation to their experimental realization (full details in ref (25)). It assumes the following: particles are single domain, in a size range of ∼10–100 nm; Brownian, not Neél, is the domain relaxation mechanism; remnant magnetization is present brought about by some variety of magnetic ordering (e.g., ferro/ferri); and this remnant magnetization can be presumed to have a fixed orientation if a high degree of magneto-crystalline anisotropy is present.

Molecular Dynamics

Extensive computer simulations were performed, facilitated by the use of the ESPReSSo 3.3.0 simulation package, to compute the zero-field magnetic response of our low-density superball fluids.[50] Our superballs are magnetic, conveniently therefore, no explicit coupling between the dispersed superballs and suspending fluid is present. As such, we have performed bulk, three-dimensional molecular dynamics simulations using a Langevin thermostat to affect constant-temperature conditions. The thermostat introduces a velocity-dependent friction term and includes the application of random jolts to the particles. These mechanisms serve as the drag force and the Brownian motion induced by the surrounding fluid. The implementation of this method in combination with our particle model is discussed in ref (25).

In this contribution, however, we are considering bulk, three-dimensional systems at ambient temperatures; a set of criteria that comes with its own distinct challenges. To this end, simulations were conducted in a cubic box with periodic boundary conditions applied, mimicking the conditions present in bulk fluid. The central concern was how to appropriately compute the dipolar interactions, as the difficulties in dealing with these long-range forces are well documented. To handle this issue, we employed the particle-particle particle-mesh (P3M) algorithm, which efficiently handles the computations required.[51]

All of our simulations were conducted in a set of dimensionless units (as defined in ref (25)), which are denoted by a superscript asterisk. For our current investigation, we set temperature to unity T* = kT = 1, which reduces the magnetic coupling parameter as follows λ = (m*)2/T* = (m*)2. The equations of motion were propagated using the velocity Verlet algorithm, where a time step of Δt* = 0.002 was found to be the most suitable choice. The accuracy of the P3M algorithm is determined by the absolute root-mean-square error in the dipolar forces. This was set as ΔF* ⩽ 10–4 for all simulations conducted. In this study we are solely concerned with equilibrium properties; the values of mass, inertia, and translational/rotation friction coefficients are physically inconsequential. Each was set to unity with the exception of the rotational friction coefficient, which takes the value of a third. This choice provides the quickest possible convergence to equilibrium.

The protocol used for each simulation consisted of the following steps. An initial configuration of N = 512 monodisperse superball particles was generated. Particle positions were randomly distributed over the volume of the simulation box, and dipole orientations were randomly distributed over the complete solid angle. The system was then propagated through an equilibration cycle. The number of steps required to achieve equilibration was typically on the order of 2–9 × 106 Δt*. The sizable range quoted is due to the strong dependence of the equilibration time on the exact physical conditions. Magnetic fluctuations, as well as the system energy, were monitored to ensure a steady state was achieved prior to any measurements being taken.

Once equilibration was achieved, a production period of 200 blocks was performed, with block-lengths ranging from 1–3.5 × 104 Δt*. A single block included measurements of the quantities of interest, e.g., magnetization and energy, evenly spaced at intervals of 100 Δt*. The ISMS was calculated using block averaging, which allowed an appropriate estimate of the associated error (standard deviation of the mean) to be made. Furthermore, on completion of a single block, a snapshot of the system configuration (particle positions and dipolar orientations) was recorded to facilitate the subsequent structural analysis. The statistical independence of the measurements was verified by careful consultation of the |M|2 autocorrelation function.