Open Access Article

This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

DOI: 10.1039/C7SM00796E
(Paper)
Soft Matter, 2018, Advance Article

Armin Kögel^{a},
Pedro A. Sánchez^{b},
Robin Maretzki^{a},
Tom Dumont^{a},
Elena S. Pyanzina^{c},
Sofia S. Kantorovich^{bc} and
Reinhard Richter*^{a}
^{a}Experimentalphysik 5, University of Bayreuth, 95440 Bayreuth, Germany. E-mail: reinhard.richter@uni-bayreuth.de
^{b}University of Vienna, Sensengasse 8, Vienna, 1090, Austria
^{c}Ural Federal University, Lenin av. 51, Ekaterinburg, 620000, Russia

Received
22nd April 2017
, Accepted 8th December 2017

First published on 11th December 2017

We investigate the phase separation of a shaken mixture of glass and magnetised steel spheres after a sudden quench of the shaker amplitude. After quenching, transient networks of steel spheres emerge in the experiment. For the developing network clusters we estimate the number of spheres in them, and the characteristic path lengths. We find that both quantities follow a log-normal distribution function. Moreover, we study the temporal evolution of the networks. In the sequence of snapshots we observe an initial regime, where the network incubates, followed by a temporal regime where network structures are elongated and broken, and finally a regime where the structures have relaxed to compact clusters of rounded shapes. This phaenomenology resembles the initial, elastic and hydrodynamic regimes observed by H. Tanaka [J. Phys.: Condens. Matter, 2000, 12, R207] during the viscoelastic phase separation for dynamically asymmetric mixtures of polymers. In order to discriminate the three regimes we investigate in the experiment order parameters like the mean number of neighbors and the efficiency of the networks. In order to capture the origin for a viscoelastic phase separation in our granular mixture, we use a simple simulation approach. Not aiming at a quantitative description of the experimental results, we rather use the simulations to define the key interactions in the experimental system. This way, we discover that along with dipolar and steric interactions, there is an effective central attraction between the magnetised spheres that is responsible for the coarsening dynamics. Our simulations show as well three regimes in the evolution of characteristic order parameters.

In the last decade, networks that are neither completely regular, nor completely random have attracted much interest.^{1} These networks are highly clustered like regular lattices, yet have small characteristic path lengths like random graphs. Examples comprise the connections in-between routers of the internet,^{1} citation networks^{2} or transportation networks in plant leaves,^{3} slime mould,^{4} or rail systems.^{5} These networks may evolve in time, but have no transient nature.

In contrast, transient networks emerge during the phase separation of dynamically asymmetric mixtures, like a solution of a polymer in a less viscous solvent.^{6} Here, the slow component is the polymer, the fast one is the solvent. For these kind of dynamically asymmetric mixtures a new type of viscoelastic phase separation has been proposed by Hajime Tanaka.^{7} This model complements the solid model and the fluid model of phase separation, which are both well known. In the viscoelastic model the phase separation starts with an initial regime – the incubation of a network – where an order parameter like the dominant wave number only slowly changes. This is followed by an elastic regime, where network like structures are elongated and eventually broken. This regime is characterized by a drastic change of the order parameter. In the final, hydrodynamic regime, “the network like structure relaxes to a structure with rounded shape and the domain shape starts to be dominated by the interface tension as in usual fluid–fluid phase separation”.^{7}

Here, we make the supposition that a dispersion of magnetic particles in a magnetically neutral granular medium may also be regarded as a dynamically asymmetric mixture. In such a system, the magnetic interactions mediate the elastic forces.

Indeed, it is well known from ferrofluids, a colloidal dispersion of magnetic nanoparticles,^{8} that they form chains and aggregates. This is especially true for a dispersion of cobalt nanoparticles.^{9} The latter excel by their high magnetic dipole moment. In the nanometer scale, however, the formation of networks is difficult to investigate because the imaging relies on transmission electron microscopy. The latter requires frozen samples^{10} which impedes to observe the dynamics. In contrast, in the micrometer range, networks of ferromagnetic particles and their dynamics can be inspected by optical means. In this way, Wen et al.^{11} observed the aggregation dynamics in a suspension of magnetised, nickel-covered glass spheres (diameter d = 50 μm) and compared the observed structures with numerical simulations based on dipole–dipole-interactions without noise. In the same range, Snezko et al.^{12} investigated the structure formation of electromagnetically driven Ni-particles. They observed compact clusters, rings, chains and rudiments of networks. Stambaugh et al.^{13,14} experimented with spheres in the cm-range, with permanent magnets of different strength incorporated in them. They investigated pattern formation and segregation. Due to the small size of the experimental systems, however, networks were not reported.

For our topic, the experiments by Blair and Kudrolli^{16} are path-breaking. The authors vibrated a mixture of permanently magnetised steel spheres and glass beads by means of a shaker. They established a phase diagram under variation of the acceleration amplitude and the filling fraction. It shows a hysteretic transition between a homogeneous gas phase and a phase, where single particles and agglomerates coexist. Quenching the system from the gas phase (cf. our snapshot in Fig. 1) into the mixed phase, the authors found loose networks like in Fig. 2a and compact networks, like in Fig. 2b.

Fig. 2 Examples for networks, as formed 5 seconds after the quench, for Γ = 1.8 g (a) and Γ = 1.9 g (b). After ref. 15. A movie of the network formation is provided as supplementary movie SM1 (ESI†). |

We revisit this rich system in order to tackle for the first time the question, whether network formation in a mixture of magnetised beads in a magnetically neutral phase can be captured by the model of viscoelastic phase separation. For our study we exploit the fact that the granular temperature can conveniently be quenched by switching the amplitude of the vibration exciter. This is a clear advantage comparing to ferrofluids, where lowering the temperature is slow and may also induce a freezing of the solvent.^{17}

The article is organised as follows. In the next section (Section 2) we describe the basics of the experiment. Thereafter we present our experimental results (Section 3), which comprise a network analysis (Section 3.1), and the temporal evolution of characteristic measures derived therefrom (Section 3.2). Then we describe our simulation method (Section 4). After that we analyse the simulated structures (Section 5) and compare them with the experimental ones. In this way, we are able to understand the nature of the magnetic interaction and the part of the glass beads. Finally (Section 6), the summary of our investigation is provided and the missing aspects are outlined.

Unavoidably, the electromagnetic vibration exciter generates a magnetic stray field. Its vertical component B_{z} has been measured by means of a Hall probe (Lakeshore, type MNA-1904-VH) and a teslameter (Lake shore, type 450). The results are shown in Fig. 4. The solid green line marks the radial dependence B_{z}(r) at a height of z = 10 cm, which has a maximum in the centre. The red dots give B_{z}(z) at r = 0 cm. The values decay drastically with increasing z, like the field of a dipole (dashed line). In order to minimize the influence of the stray field onto the steel spheres, the vessel is fixed via a hollow brass rod (outer diameter 20 mm) and a bearing to the flange of the exciter. In this way, the bottom of the experimental vessel is situated 53.5 cm above the top of the exciter, where B_{z} < 61 μT. The latter is in the range of B_{z} of the earth in central europe (44 μT) and shall be neglected.

The vessel is filled with a mixture of steel spheres and glass beads. Their position is recorded by means of a charge-coupled-device camera (Lumenera Corporation, type LU135M) with 8-bit resolution, where each frame has 1392 × 1040 pixel. The spheres are illuminated from below, by an electroluminescent film (Zigan Displays) mounted beneath a glass plate, as sketched in Fig. 5. The dimensions (width × length × height) of the rectangular vessel are 200 × 285 × 12.85 mm^{3}. The vessel has no lid. In order to achieve a light-weight but stiff vessel its bottom plate is made from aluminum honeycomb material (CEL Components S.R.L.). During the experiments we record the acceleration amplitude Γ with an acceleration sensor (Brüel & Kjaer, type 4509 B002), which is mounted beneath the vessel.

H_{i} = H − DM.
| (1) |

M(H) = ±M_{R} + χ_{0}H_{i},
| (2) |

Fig. 6 Magnetisation M vs. internal magnetic field H_{i} following eqn (1) of the selected steel spheres. The red rectangle in (a) marks the range of the zoom displayed in chart (b). Here the filled circles give the measured data, the solid lines present fits by eqn (2) for estimating the magnetic remanent magnetisation M_{R} and the susceptibility χ_{0}. |

Fig. 7 Remanent dipole moment μ_{R} (l.h.s.) and coercivity H_{C} vs. applied magnetic induction for the selected steel spheres. |

Material | Soda-lime glass | Steel 1.3505 |
---|---|---|

Manufacturer | Sigmund-Lindner.com | Isometall.de |

Type | Type P | DIN5401G10 |

Radius r (mm) | 2.0 ± 0.02 | 1.5 ± 0.02 |

Volume V (mm^{3}) |
33.51 | 14.14 |

Density ρ (kg m^{−3}) |
2580 | 7610 |

Mass m (g) | 0.084 | 0.108 |

Remanence M_{R} (kA m^{−1}) |
— | 340 ± 3 |

Susceptibility χ_{0} |
— | 76 ± 2 |

Coercivity H_{C} (kA m^{−1}) |
— | 4.46 ± 0.3 |

Dipole moment μ_{R} |
||

(10^{−4} A m^{2}) |
— | 4.83 ± 0.02 |

The filling fraction of the glass spheres, ϕ_{g}, as well as the one of the magnetised spheres, ϕ_{m}, can be varied in the experiment. Following ref. 16 we define the filling fraction via

(3) |

Two neighbouring magnetised steel spheres attach to each other if they are close enough and properly oriented. After forming a dimer, the mass of the new cluster has doubled compared to each particle, which lowers its mobility. Moreover, its magnetic moment has increased as well. Consecutively, further particles tend to be attracted to the ends of the dimer, creating a chain of three magnetised spheres, whereas a triangular configuration is hardly observed, because it is energetically less favourable. The trimers extend to longer chains. With their length increasing, it becomes more likely that some steel spheres are also connected to the side of a chain, forming various branched structures^{22,23} as shown in Fig. 8a. With increasing length, the chains bend more easily and collisions with glass spheres become more frequent. This can lead to the formation of rings, like displayed in Fig. 8b, or to the disruption of chains. The further evolution of the pattern is very sensitive to the value of the shaker amplitude. For lower amplitude, loose networks are favoured (cf. Fig. 2a), whereas for slightly higher amplitudes, compact networks and areas with dense hexagonal packing (“crystallites”) form earlier (see Fig. 2b).

Fig. 8 Clusters which emerge can reduce their magnetic energy by forming branched chains (a) and rings (b). |

In the following Section 3.1 we introduce and apply different network measures to the ferrogranular structures. These measures are then exploited in Section 3.2 to describe the temporal evolution of the transient networks.

3.1.1 Number of neighbours. A basic measure for networks is the number of neighbours k to which each node is connected by edges. It is in a two-dimensional layer geometrically limited to 6. Next we exemplary analyse the two states illustrated in Fig. 2 by the distribution of k, as shown in Fig. 9. For Γ = 1.8 g (bright), almost all spheres are connected to neighbours, whereas at Γ = 1.9 g (dark), some spheres remain alone (k = 0). For the lower amplitude, the numbers of spheres which are situated within a chain or a ring (k = 2) or at chain ends (k = 1) is about twice as high as those for a slightly higher amplitude. This is in agreement with Fig. 2b where crystallites are more frequent than in (a). The greater compactness of the state of Fig. 2b is also reflected in the mean number of neighbours^{1} , as presented in Table 2.

Fig. 9 Relative occurrence of the number of neighbours (k) for Γ = 1.8 g (bright) and 1.9 g (dark) and ϕ_{g} = 0.15, ϕ_{m} = 0.15. |

Acceleration Γ (g) | |
---|---|

1.8 | 2.39 ± 0.06 |

1.9 | 3.15 ± 0.10 |

3.1.2 Number of spheres in the clusters. The cluster size at a specific time depends strongly on the acceleration amplitude Γ. As seen in Fig. 2, larger clusters are favoured by a slightly higher Γ, that enhances the mobility of the still solitary spheres.

where μ denotes the mean value and σ the standard deviation. It describes the normal distribution of logarithmically distributed random numbers and is ubiquitous in nature^{24} for systems with positive definite observables. This function was successfully applied to the particle size of smashed objects^{25} and may fit as well for repeatedly smashed clusters. Moreover it is also known from random network structures.^{26}

where N_{T} denotes the total number of clusters, and erf the Gauß error function. The mean m and the standard deviation s are related to the log-normal distribution via^{24}

The border of the clusters, as marked by the red lines in Fig. 10, are automatically detected. The number of spheres within the lines yields the cluster size. Many individual spheres and a few large clusters are dominating the picture, where the size of the latter is limited by the finite size of the vessel. Fig. 11 presents the distribution of the cluster sizes for Γ = 1.8 g five seconds after the quench. The distribution resembles the log-normal distribution

(4) |

Fig. 10 Detected clusters in the network depicted in Fig. 2b. Connected steel spheres are bordered by the red line. The inner voids, created by loops, are marked by the blue line. |

Fig. 11 Distribution of cluster size for Γ = 1.8 g, and ϕ_{g} = 0.27, ϕ_{m} = 0.25, determined from 50 frames. |

For a careful comparison of experiment and eqn (4), we plot in Fig. 12 the cumulative distribution of the cluster size for two different accelerations. The solid lines display fits of the distribution function

(5) |

(6) |

(7) |

Fig. 12 Distribution function fitted to the cumulative incidence of the cluster size for Γ = 1.8 g and 1.9 g, with ϕ_{g} = 0.27, ϕ_{m} = 0.25. |

In Fig. 12 the measured data for Γ = 1.8 g can well be fitted by eqn (5), whereas at Γ = 1.9 g deviations are prominent at small cluster size, which is reflected in Table 3.

Acceleration Γ (g) | 1.8 | 1.9 |
---|---|---|

m_{real} |
16.0 | 14.2 |

m_{fit} |
16.5 | 185 |

s_{real} |
27.0 | 70.6 |

s_{fit} |
37.5 | 550 |

3.1.3 Characteristic path length. Fig. 13 presents a sight of Fig. 2b, where all spheres were replaced by nodes, and the edges connecting the spheres by lines. From such a reduced graph the characteristic path length

can be defined, where N counts the number of nodes or spheres, and d_{ij} measures the shortest path length in-between the i-th and the j-th node. It measures the minimum number of edges one has to trespass in order to go from node i to node j.^{1} Therefore, the characteristic path length L is an average of all shortest path lengths within a cluster.

(8) |

Fig. 13 Reduced graph of Fig. 2b, where only the connecting lines in-between the steel spheres are shown. The supplementary movie SM2 (ESI†) is showing the evolution of the reduced network graph. |

In order to determine L of every cluster, we use the software package networkX^{19} for the script language python.^{20} Similar to the previous section, we investigate again the cumulative distribution of the cluster size, now measured in L. As shown in Fig. 14, the cumulative distributions for both acceleration amplitudes are well fitted by the integral of the log-normal-distribution eqn (5). Comparing the mean of the measured (real) and the fitted (fit) distributions in Table 4 shows only a tiny difference, which confirms the quality of the fit.

Fig. 14 Cumulative incidence of clusters vs. the characteristic path length for Γ = 1.8 g and 1.9 g, with ϕ_{g} = 0.27, ϕ_{m} = 0.16. The solid lines mark fits by eqn (5). For fitting parameters see Table 4. |

Acceleration Γ (g) | 1.8 | 1.9 |
---|---|---|

m_{real} |
4.55 | 7.80 |

m_{fit} |
4.58 | 8.46 |

s_{real} |
5.25 | 12.1 |

s_{fit} |
6.39 | 23.1 |

So far we have characterized the emerging network by means of the number of neighbours k and by the characteristic path length L. We have focused on these quantities because they serve in the following to characterize the evolution of the network. However, its worth to note, that the distribution of clusters measured by means of the gyration radius shows as well a log-normal distribution. Moreover, also the distribution of the loop sizes is an interesting observable. An analysis of the latter observables will be given in a subsequent publication.

3.2.1 Evolution of the mean number of neighbours. The mean number of neighbours is plotted in Fig. 15 vs. the time for two different acceleration amplitudes. We start in the gas phase where = 0, and quench at t = 0 s. During the first five seconds after the quench increases rapidly. During this time the steel spheres form first dimers and trimers. During the further evolution grows more slowly.

which describes many limited growth processes in nature very well.^{27} In eqn (9) _{max} denotes the saturation value of the growth and t_{0} the time when equals to the half of the saturation value. With increasing t_{0} the network grows more slowly. The exponent p describes the curvature of the graph. For p ≫ 1 the logistic function increases in the beginning more slowly, but grows more rapidly around t_{0}, where it approaches fast its saturation value.

Fig. 15 Temporal evolution of the average number of neighbours for two different accelerations Γ, and ϕ_{g} = 0.27, ϕ_{m} = 0.16. The solid lines indicate fits by eqn (9). |

In a first attempt, we try to characterise the temporal evolution of by the logistic function

(9) |

The fit of (t) agrees rather well with the experimental values, as displayed in Fig. 15. Table 5 presents the parameters derived from the fit. For higher Γ the saturation value _{max} increases. This is in agreement with our qualitative observations, which relate a higher Γ with more compact structures. Also t_{0} varies strongly with Γ. For Γ = 2.0 g half of the saturation value is reached at t_{0} = 1.9 s, whereas for Γ = 2.1 g this happens at t_{0} = 3.4 s, which is considerably later.

Acceleration Γ (g) | 2.0 | 2.1 |
---|---|---|

_{max} |
3.3 | 3.8 |

t_{0} (s) |
1.9 | 3.4 |

p | 0.75 | 0.86 |

3.2.2 Evolution of the efficiency. For the characteristic path length L (8) the length in-between spheres on separate clusters diverges. Therefore L must be estimated for the individual clusters, and later be averaged. This can be circumvented by using the efficiency

of a network.^{28} Here, N counts the number of spheres, and d_{ij} is the shortest path length in-between sphere i and sphere j. In this way E measures the parallel information exchange in-between the nodes of the network.^{28}

(10) |

The efficiency E is plotted in Fig. 16 vs. time. For all investigated values of Γ, we observe a monotonous increase of E. In the first five seconds after the quench at t = 0 s E increases rapidly, which is followed by a more moderate evolution afterwards. Here, small steps in the evolution occur if two larger clusters are merging. This is an outcome of the finite size of the container and the limited statistics. In contrast to the evolution of the average number of neighbours, (see Fig. 15 above) the curves in Fig. 16 do not indicate a saturation. Therefore, a fit by the logistic function is not convincing. We note in Fig. 16 as well that higher values of E are related with a higher acceleration amplitude Γ.

Fig. 16 The efficiency E according to eqn (10) vs. time t for three different acceleration amplitudes Γ. For all Γ the filling fractions (3) are ϕ_{g} = 0.15, ϕ_{m} = 0.20. |

It is interesting to know how the filling fraction ϕ will influence E(t). In Fig. 17, we plot E(t) for four different values of ϕ_{m}. For each Γ we observe a monotonous increase of E, which is, similarly to that in Fig. 16, first steep and then moderate. For higher values of ϕ_{m}, one may expect larger clusters and larger path length L. Due to the limit of k = 6 this will result in a smaller efficiency. Indeed, at e.g. t = 50 s, a higher value of ϕ_{m} leads to a lower E for all filling fractions. However at t ≈ 120 s we observe a crossover of the curves for ϕ_{m} = 20% and 25%. This may be due to the limited statistics of our experimental arrangement.

Fig. 17 The efficiency E according to eqn (10) vs. time t for Γ = 1.9 g. The filling fraction for the glass beads was fixed to ϕ_{g} = 0.15, whereas the one for the magnetised beads was varied according to ϕ_{m} = 0.15, 0.20, 0.25, and 0.30. |

The graphs of (t) (cf. Fig. 15) as well as E(t) (see Fig. 16 and 17) clearly display two distinct time scales. A short time scale ≈2t_{0} of a few seconds, where the network is emerging, and a long time scale ≈100·t_{0} where its structure is being transformed. These two time scales indicate an initial and an elastic regime, predicted in ref. 7 for viscoelastic phase separation. Moreover, the emergence of compact hexagonal arrangements of magnetised spheres can clearly be discriminated in the supplementary movies SM1 and SM2 (ESI†). However, they are not yet dominating the pattern, as it is predicted for the hydrodynamic regime.^{7} Consequently, it is difficult to unveil the transition in the evolution of the order parameters in the Fig. 15 and 16.

For times larger than 180 s we can detect an increase in E. However, we do not show this here, because we rely not fully on the statistics due to finite size effects, i.e. collisions of the clusters with the edge of the vessel. The latter are becoming increasingly important at long times, where the cluster size is in the range of the experimental vessel. Note that already 5 s after the quench networks can approximate the diameter of the vessel, as can be seen in Fig. 2b.

Due to mechanical limitations it is not straight forward to enlarge our experimental setup by an order of magnitude. Instead we re-examine the long term behaviour in a simple numerical model.

Before introducing our theoretical approach, it is important to underline that it does not aim at any quantitative characterisation of the system. Instead we focus strictly on establishing the basic ingredients that can reproduce qualitatively the main properties of the network evolution, particularly the ones corresponding to slow dynamics.

A simple starting point to define a minimal model of spheres that exhibit magnetic interactions is to consider that they carry a magnetic point dipole in their centres. Despite this description is only exact for homogeneous monodomain spheres, it has been used extensively as a useful approximation for the study of physical systems with very different characteristic sizes, ranging from magnetic nanoparticles to systems of macroscopic spherical magnets. The highly anisotropic nature of the dipole–dipole interaction leads to a strongly directional spontaneous self-assembly of the spheres, that tend to form linear chains with head-to-tail arrangements of their dipoles. This effect is well known in the study of ferrofluids and other systems of ferromagnetic nanoparticles.^{29–34} However, the networks of magnetised steel spheres observed in our experiments include a significant fraction of close packing crystallites, a structure that is very unfavourable for a system driven only by dipole–dipole self-assembly.

One can hypothesise that the dipolar approximation is still valid for the magnetic interactions in this system and that the close packing topology has its origin in the kinetic effect of the bath of non magnetic beads. Alternatively, one may claim that the size of the steel spheres and their experimental magnetisation curve described by eqn (2) suggest that they may actually have a multidomain nature, making the dipolar approximation too inaccurate. In this case, the magnetic dipole moment of each individual sphere would depend in part on the field created by its surrounding neighbours. The overall effect of this complex interplay would be to make the effective interaction between the steel spheres less anisotropic than the corresponding to pure dipolar spheres, allowing the formation of close packing crystallites. In order to ellucidate which is the correct hypothesis, we chose a computer model that allows us to easily separate such effects, while keeping a minimalistic approach. This model and the simulation method are described in the following sections.

(11) |

(12) |

Our strategy to clarify what is the actual drive of the coarsening dynamics observed in our experiments is to perform simulations for three different systems, corresponding to distinct qualitative combinations of the interactions defined above, and to compare the resulting networks to the experimental ones. Following an order of increasing complexity, first we simulated a monodisperse system of ideal soft core dipolar spheres (S1), i.e., a system of identical beads with characteristic diameter σ_{m} and magnetic moment . These beads interact via the dipole–dipole potential (12) and a purely repulsive WCA potential, obtained from eqn (11) by taking r_{cut} = 2^{1/6}σ, σ = σ_{m} and ε = ε_{r}. For the second system (S2), we add to S1 a fraction of identical non magnetic beads, representing the glass spheres, with characteristic diameter σ_{g}. The repulsion between these non magnetic beads is calculated in the same way by using σ = σ_{g}, whereas for the repulsion between the magnetic and non magnetic beads we take σ = (σ_{m} + σ_{g})/2. Finally, the last system (S3) is obtained by adding to the interactions defined for S2 the aforementioned isotropic attraction between the magnetic beads, also obtained from eqn (11) by taking r_{cut} = 2.5σ, ε = ε_{a} and σ = σ_{m}.

m_{i}(d_{i}/dt) = _{i} − Γ_{T}_{i} + _{i,T}
| (13) |

_{i}·(d_{i}/dt) = _{i} − Γ_{R}_{i} + _{i,R},
| (14) |

As it is frequently done in coarse-grained simulations, we take an arbitrary system of reduced units for all physical parameters, with the only constraint of keeping the relative scales of the relevant quantities as close as possible to the experimental ones. Mass units are defined by the mass of the steel and glass beads, m_{m} and m_{g} respectively. These were chosen to be equal, so that we take m_{m} = m_{g} = 1. Lengths are measured in terms of the respective diameters of the beads. Keeping the experimental ratio, we take σ_{m} = 3 and σ_{g} = 4. Energies are measured in the scale of the soft-core repulsive interactions, so that we chose ε_{r} = 10. Finally, in general one can define formally the reduced time scale in the system in terms of the units of mass, length and energy. However, here one has to take into account that the temperature plays a role that is only defined qualitatively by comparison with the experimental results. This makes it more accurate to also determine the time scale of the simulations by fitting the coarsening dynamics of the simulated systems to the experimental measurements. It has also the additional advantage of letting the choice of the rest of dynamic parameters to be arbitrary. Therefore, we can choose the reduced values of the friction constants to be Γ_{T} = 1 and Γ_{R} = 3/4, which are values known for providing a fast relaxation in this type of simulations.^{37,38} Finally, in order to ensure isotropic rotations, we set the inertia tensors, _{i}, as the identity matrix.

In order to mimic the experimental setup as close as possible without suffering its finite size effects we simulate a pseudo-infinite two-dimensional system by fixing the z coordinate of all particles to zero and using lateral periodic boundaries for the x–y plane. Under these conditions, the calculation of the long-range dipole–dipole interactions requires an efficient approximate method. For this we used the dipolar-P^{3}M algorithm,^{39} combined with the dipolar Layer Correction method^{40} in order to take into account the slab geometry of our system. The simulations were carried out with the package ESPResSo 3.2.0.^{41,42}

For the sake of simplicity, in simulations we focus on a single set of experimental parameters to fit the model and analyse the coarsening dynamics. Specifically, in all the simulations the area fraction of the magnetic beads is ϕ_{m} ≈ 0.18, whereas in systems S2 and S3 the area fraction of non magnetic beads is ϕ_{g} ≈ 0.15, so that ϕ_{m} + ϕ_{g} ≈ 0.33. These values were obtained by placing 1466 magnetic and, for S2 and S3, additionally 687 non magnetic beads in a square simulation box with side length L = 240. In all cases the temperature during the first simulation step was taken to be T_{1} = 2. Finally, in order to fit the simulation model to the experimental results, we explore diverse values of the quenching temperature, T_{2} = {0.3,0.5,0.75}, energy scale of the isotropic attraction in S3, ε_{a} = {0.3,0.5,1}, and dipole moment of the magnetic beads, μ = ||. The values of μ were chosen to provide reasonable values of the so-called dipolar coupling parameter, λ, at the sampled values of T_{1} and T_{2}. This parameter, frequently used to characterise the self-assembly behaviour in systems of magnetic dipolar particles, is defined as the ratio between the minimum dipole–dipole energy of a pair of particles and the thermal energy of the system. Here, we can express it as:

(15) |

According to this definition, the sampled values of μ correspond to the intervals λ ∈ [0.5,1.75] for T_{1} and λ ∈ [2,7] for T_{2}. However, in the following we will focus the discussion of the results only on the most interesting values of these parameters. As a last comment on the details of the simulation protocol, at least 16 independent simulation runs were performed for every set of parameters to obtain reasonably good statistics. These simulation runs consisted of 3 × 10^{5} integration steps at T_{1}, with a simulation time step δτ_{1} = 0.0005, and the same amount of integration steps at T_{2}, with a simulation time step δτ_{2} = 0.01.

In order to perform a quantitative characterisation of the simulated networks, we first have to identify the clusters of connected magnetic particles. For this we use a simple distance criterium: two magnetic particles are connected if their centre-to-centre distance is less than 1.3σ_{m}. Note that this threshold corresponds to applying the same factor used for identifying the experimental networks, 1.13, to the distance at which the soft core steric repulsion used in the simulations vanishes, r_{cut} = 2^{1/6}σ_{m}. Once the networks are identified, their parameters are calculated in exactly the same way as we did for the experimental data.

Fig. 18 Snapshots of the simulations of mixtures of magnetic and non magnetic particles obtained for λ = 5 and T_{2} = 0.5, at different integration steps n_{δτ}, after quenching: n_{δτ} = 2500 (a0–a3), n_{δτ} = 10^{5} (b0–b3) and n_{δτ} = 2.5 × 10^{5} (c0–c3). The column (a0, b0, c0) corresponds to the system without central attraction, S2 (ε_{a} = 0), and the rest to different strengths of the central attraction, S3: ε_{a} = 0.3 (a1, b1, c1), ε_{a} = 0.5 (a2, b2, c2), and ε_{a} = 1.0 (a3, b3, c3). Magnetic particles are plotted in black, non magnetic in light grey. Neighbouring contacts between particles forming clusters are marked in orange (see the main text for the criterium to identify neighbouring particles). Supplementary movies (ESI†) are provided for S2 (SM3) as well as for S3, ε_{a} = 0.5 (SM4) and S3, ε_{a} = 1.0 (SM5). |

For systems S2, as shown in the examples of column (a0, b0, c0), one can observe that the magnetic beads form structures composed by short linear chains, rings and branched chains. These structures evolve with time, albeit without obvious qualitative changes: Only the size of clusters grows slightly and some of them merge to form side-by-side arrangements. Still the general structure is much more loose than the one observed in the experiments.

In contrast, introducing a central attraction makes a strong impact on the overall structure: even for a weak central attraction—ε_{a} = 0.3, shown in column (a1, b1, c1)—one sees larger clusters with bundles of chains and/or rings that deplete/enclose the non magnetic particles. The growth of these clusters with time is also clearly stronger. This trend becomes even more pronounced as the central attraction grows, as is shown in the columns (a2, b2, c2) for ε_{a} = 0.5 and (a3, b3, c3) for ε_{a} = 1. For the latter case, we can see how the system separates very early into magnetic and non magnetic particles, with the magnetic ones forming compact clusters of nearly hexatic packing. This observation already suggests that the central attraction may be an essential ingredient to capture the dynamics of network formation in experiments. A more rigorous comparison between simulation and experimental results is presented in the next section.

Fig. 19 shows the best fit to the experimental data (marked by □) of the values of obtained from simulations of five systems with different interactions. These are systems S1, S2 and S3 with three different strengths of central attraction, ε_{a} = {0.3,0.5,1.0}. In all cases, time t = 0 corresponds to the moment at which the quenching takes place. The fit consisted of a rescaling of the simulation time to match the experimental evolution profile obtained during the first 20 seconds after quenching. Specifically, we obtained the relationship t = 0.438τ, where t is the physical time and τ is the simulation time. This time rescaling will be applied to all remaining simulation results.

The comparison of all the simulation curves with the experimental data clearly indicates that systems S3 are much closer to the experiment than S1 and S2, which underestimate significantly. Among the S3 systems, the one with ε_{a} = 0.5 shows an excellent agreement with the experimental results within the range used for the fitting, whereas lower and higher strengths of this interaction tend to underestimate and overestimate , respectively. One can also see that, even for the best set of model interactions, systematic deviations from the experimental values tend to increase with time. This can be explained by the growing impact of finite size effects on the experimental measurements. Therefore, we conclude that the presence of the central attraction is an essential requirement to capture the experimental dynamics in our minimalistic modelling approach and, in particular, that the intermediate value ε_{a} = 0.5 provides the best approximation to the experimental results obtained for shaking amplitude Γ = 1.8 g. For simplicity, from now on we will refer to the model system S3 with ε_{a} = 0.5, λ = 5 and T_{2} = 0.5 as S3*.

Once we determined the best model for the experimental system of reference, we can examine in simulations the long time properties of the networks of magnetic beads under ideal measurement conditions. In the next sections we present the simulation results for the average number of neighbours and the network efficiency.

The existence of three dynamic regimes, as evidenced by the behaviour of the parameter , suggests that the system may indeed display a viscoelastic phase separation, so that such regimes may correspond to the ones proposed by Tanaka:^{7,43} initial, elastic and hydrodynamic, respectively. However, their accurate interpretation requires subtle considerations. As pointed out above, regime (i) corresponds to the aggregation of isolated particles into clusters whose structure is dominated by the dipolar term of the pair interaction. This leads to networks with mainly chain-like topologies. The transition from regime (i) to regime (ii) signals the point where the fraction of particles with more than two neighbours, favoured by the central part of the pair interaction, starts to be significant. This suggests that the main phenomena in region (ii) may correspond to changes in the topology of the networks, that evolves from chain-like to hexatic crystallites. Finally, the interpretation of the transition between regimes (ii) and (iii) is hard to establish without further elements of analysis. These are provided in the next section, where results for the network efficiency and the mean cluster size are discussed.

A complementary picture of the different dynamic phenomena present in the system can be obtained by computing a third parameter: the mean relative cluster size, . This is simply the average number of magnetic beads per cluster divided by the total number of magnetic beads in the system. Fig. 21(b) shows the simulation results for this parameter, together with the two sets of regime limits determined for (indicated by the borders of the grey area) and E (marked by dashed lines). We can see that initially there is a very slow growth of , that remains up to an intermediate point between the beginning of regime (ii) of the network efficiency and the end of regime (ii) for the number of neighbours. After such point, the growth rate of becomes much larger. This fast growth is a clear signature of the merging of separate clusters into larger ones.

Finally, the comparison of the results shown in Fig. 20 and 21 allows us to complete the sequence of main dynamic phenomena governing the evolution of the networks. After the initial interval in which a fast aggregation of the magnetic beads into clusters with a mainly chain-like structure takes place, driven by the long range dipolar interactions and corresponding to regime (i) of the evolution of the mean degree, the formed clusters experience a compaction by slowly replacing their chain-like morphologies by hexatic crystallites allowed by the short range central attraction. This compaction is favoured by the collisions of the non magnetic beads on the clusters. We identify this regime with the elastic phase of Tanaka's viscoelastic transition. Finally, the slowest phenomenon is the merging of separate clusters into larger ones, which becomes the dominant effect near the end of the regime (ii) of and manifests itself in the fast growth of both, and E. We consider this to correspond to the hydrodynamic regime of Tanaka's transition. The fact that the merging of the clusters, producing the final phase separation of magnetic and non magnetic beads, happens mainly at very long times is the consequence of the slow mobility of the clusters compared to individual particles: recently, it has been shown that such mobility tends to decrease with the cluster size.^{44} As a final observation, the slight slowing down of the growth of E corresponding to its regime (iii) can be explained as the onset of a saturation of the impact of the merging process in this parameter: whenever two separate clusters merge, all the pairs of beads belonging to each one of the original clusters start to contribute to increase E. However, according to the definition (10), when the merged clusters are large most of the newly networked pairs have a little contribution to this growth because their shortest path lengths are anyway very large. In other words, statistically, the merging of many small clusters into medium sized ones has a stronger impact on E than the merging of few large clusters into very large ones.

Fig. 21 Long time evolution of the network efficiency, E, obtained for simulations with model system S3* (symbols), and semilogarithmic fits to the three regimes observed, labelled as (i) to (iii) (solid lines). Dotted vertical lines indicate the limits between regimes. The slopes of the fits are: 0.14 (i), 0.47 (ii) and 0.31 (iii). (b) Long time behaviour of the mean relative cluster size, , for the same system. To ease the comparison with Fig. 20 and 21(a), the regime limits shown in such figures are also indicated here by the dark background and the vertical dotted lines. |

Still it is safe to conclude that experiment as well as simulations show an initial regime (i), where short chains of dipoles are emerging, which is followed up by a time span exhibiting the formation of a complex and contracting network of chains, and finally (iii) a regime characterised by compact clusters of magnetised spheres with rounded shape, which slowly evolve. We relate these states with the initial (i), elastic (ii) and hydrodynamic (iii) regimes, first observed by Tanaka^{7} during the phase separation of dynamically asymmetric mixtures of polymers. In our granulate, the magnetised beads represent the slow component, whereas the glass beads act as the fast one.

The transient networks could be characterised in the experiment by the distribution of cluster size, as measured in number of spheres, and the mean average path length. Interestingly, at short and moderate time-scales, these distributions can well be fitted by log-normal functions.

In search for a suitable order parameter of the viscoelastic phase separation, we studied the mean degree of a node, the efficiency of the network and the average cluster size. In experiment and simulation the short time evolution of the degree is well described by a logistic function, that is characteristic for limited growth processes. In the simulations, we can access also long-time behaviour of the mean node degree, that clearly shows three distinct regimes, thus, making this observable a good candidate to serve as an order parameter for our system. Moreover, the growth of the efficiency and the cluster size exhibits different time-scales, being the former significantly larger than the corresponding to the degree. Here, the statistics needs to be improved in follow up experiments and simulations in order to better discriminate the elastic and the hydrodynamic regimes.

The comparison between the experiment and the model unveils that besides magnetic dipoles, an effective central attraction is essential for the transient nature of the network. This central attraction is a consequence of the magnetisation curve of the steel spheres. Each sphere represents a dipole with susceptibility, that we can name in short: a DipSus. In reality, the attraction in the experiment is the result of the Kelvin force, and is intrinsically anisotropic. For simplicity, our model incorporates just an isotropic approximation to this attraction. In other words, we model our DipSus system as a ferrogranular analogue of a Stockmayer fluid.^{45}

To conclude, shaken ferrogranulate exhibits interesting features, including network formation and viscoelastic phase separation. More in depth studies are necessary to find commonalities and differences to networks recently observed in systems of externally induced dipoles.^{46–48}

Last, but not least, we add a far-ranging outlook. Due to the abundance of iron in many asteroids and planets, and its high Curie temperature, the observed coarsening dynamics of susceptible and partly magnetised iron particles may play an important role in the early stage of planet formation. This process may be capable to overcome the bouncing barrier^{49} in particle agglomeration. The emerging iron “cores” may than trigger the gravitational agglomeration of the remaining matter. Related experiments under microgravity conditions are yet to be performed.

- S. Boccaletti, V. Latora, Y. Moreno, M. Chavez and D.-U. Hwang, Complex networks: Structure and dynamics, Phys. Rep., 2006, 424, 175–308 CrossRef.
- J. E. Hirsch, An index to quantify an individual's scientific research output, Proc. Natl. Acad. Sci. U. S. A., 2005, 102(46), 16569–16572 CrossRef CAS PubMed.
- E. Katifori, G. J. Szöllősi and M. O. Magnasco, Damage and fluctuations induce loops in optimal transport networks, Phys. Rev. Lett., 2010, 104, 048704 CrossRef PubMed.
- W. Baumgarten and M. J. B. Hauser, Functional organization of the vascular network of physarum polycephalum, Phys. Biol., 2013, 10(2), 026003 CrossRef PubMed.
- A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi and T. Nakagaki, Rules for biologically inspired adaptive network design, Science, 2010, 327(5964), 439–442 CrossRef CAS PubMed.
- H. Tanaka, Unusual phase separation in a polymer solution caused by asymmetric molecular dynamics, Phys. Rev. Lett., 1993, 71(19), 3158–3161 CrossRef CAS PubMed.
- H. Tanaka, Viscoelastic phase separation, J. Phys.: Condens. Matter, 2000, 12, R207 CrossRef CAS.
- R. E. Rosensweig, Ferrohydrodynamics, Cambridge University Press, Cambridge, New York, Melbourne, 1985 Search PubMed.
- Colloidal Magnetic Fluids: Basics, Development and Applications of Ferrofluids, ed. S. Odenbach, Lect. Notes Phys., Springer, Berlin, Heidelberg, New York, 2009, vol. 763 Search PubMed.
- K. Butter, P. H. H. Bomans, P. M. Frederik, G. J. Vroege and A. P. Philipse, Direct observation of dipolar chains in iron ferrofluids by cryogenic electron microscopy, Nat. Mater., 2003, 2, 88–91 CrossRef CAS PubMed.
- W. Wen, F. Kun, S. Pal, D. W. Zheng and K. N. Tu, Aggregation kinetics and stability of structures formed by magnetic microspheres, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, R4758 CrossRef CAS.
- A. Snezhko, I. S. Aranson and W. K. Kwok, Structure formation in electromagnetically driven granular media, Phys. Rev. Lett., 2005, 94, 108002 CrossRef CAS PubMed.
- J. Stambaugh, D. P. Lathrop, E. Ott and W. Losert, Pattern formation in a monolayer of magnetic spheres, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 026207 CrossRef PubMed.
- J. Stambaugh, Z. Smith, E. Ott and W. Losert, Segregation in a monolayer of magnetic spheres, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 031304 CrossRef PubMed.
- A. Kögel, Messungen zur Vergröberungsdynamik in transienten ferrogranularen Netzwerken, Bachelor thesis, Experimentalphysik 5, University of Bayreuth, 95440 Bayreuth, Germany, 2015 Search PubMed.
- D. L. Blair and A. Kudrolli, Clustering transitions in vibrofluidized magnetized granular materials, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 021302 CrossRef PubMed.
- A. V. Lebedev, Low-temperature magnetic fluid stabilized with mixed fatty acids, Colloid J., 2010, 72, 810 CrossRef.
- Open Computer Vision (openCV), http://https://github.com/opencv/opencv/wiki, accessed at 14th of September 2016.
- NetworkX: High-productivity software for complex networks, http://https://networkx.github.io/, link accessed at 12.4.2017.
- Python software foundation, http://https://www.python.org/, accessed at the 13th of September 2016, 2016.
- A. Kögel, Video of the formation of a ferrogranular network, http://www.ep5. uni-bayreuth.de/de/research/Magnetic-Soft-Matter/video/ferronetwork.html, 2015.
- T. Tlusty and S. A. Safran, Defect induced phase separation in dipolar fluids, Science, 2000, 290, 1328 CrossRef CAS PubMed.
- S. S. Kantorovich, A. O. Ivanov, L. Rovigatti, J. M. Tavares and F. Sciortino, Temperature-induced structural transitions in self-assembling magnetic nanocolloids, Phys. Chem. Chem. Phys., 2015, 17, 16601–16608 RSC.
- E. Limpert, W. Staehl and M. Abbt, Log-normal distributions across the sciences: Keys and clues, BioScience, 2001, 51(5), 341–352 CrossRef.
- T. Ishii and M. Matsushita, Fragmentation of long thin glass rods, J. Phys. Soc. Jpn., 1992, 61(10), 3474–3477 CrossRef CAS.
- J. Shackelford, The lognormal distribution in the random network structure, J. Non-Cryst. Solids, 1981, 44, 379–382 CrossRef CAS.
- D. W. Thompson, On Growth and Form, Cambridge University Press, Cambridge, UK, 2nd edn, 1942 Search PubMed.
- V. Latora and M. Marchiori, Efficient behavior of small-world networks, Phys. Rev. Lett., 2001, 87(19), 198701 CrossRef CAS PubMed.
- J. J. Weis and D. Levesque, Chain formation in low density dipolar hard shperes: a monte carlo study, Phys. Rev. Lett., 1993, 71(17), 2729–2732 CrossRef CAS PubMed.
- R. P. Sear, Low-density fluid phase of dipolar hard spheres, Phys. Rev. Lett., 1996, 76, 2310–2313 CrossRef CAS PubMed.
- R. van Roij, Theory of chain association versus liquid condensation, Phys. Rev. Lett., 1996, 76(18), 3348–3351 CrossRef CAS PubMed.
- Y. Levin, What happened to the gas–liquid transition in the system of dipolar hard spheres?, Phys. Rev. Lett., 1999, 83, 1159–1162 CrossRef CAS.
- P. J. Camp, J. C. Shelley and G. N. Patey, Isotropic fluid phases of dipolar hard spheres, Phys. Rev. Lett., 2000, 84(1), 115–118 CrossRef CAS PubMed.
- P. J. Camp and G. N. Patey, Structure and scattering in colloidal ferrofluids, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 62(4), 5403–5408 CrossRef CAS.
- J. D. Weeks, D. Chandler and H. C. Andersen, Role of repulsive forces in determining the equilibrium structure of simple liquids, J. Chem. Phys., 1971, 54, 5237 CrossRef CAS.
- D. Frenkel and B. Smit, Understanding molecular simulation, Academic Press, 2002 Search PubMed.
- J. J. Cerdà, S. Kantorovich and C. Holm, Aggregate formation in ferrofluid monolayers: simulations and theory, J. Phys.: Condens. Matter, 2008, 20, 204125 CrossRef PubMed.
- S. Kantorovich, J. J. Cerdà and C. Holm, Microstructure analyisis of monodisperse ferrofluid monolayers: theory and simulation, Phys. Chem. Chem. Phys., 2008, 10(14), 1883–1895 RSC.
- J. J. Cerdà, V. Ballenegger, O. Lenz and C. Holm, P3m algorithm for dipolar interactions, J. Chem. Phys., 2008, 129, 234104 CrossRef PubMed.
- A. Bródka, Ewald summation method with electrostatic layer correction for interactions of point dipoles in slab geometry, Chem. Phys. Lett., 2004, 400(1–3), 62–67 CrossRef.
- H. J. Limbach, A. Arnold, B. A. Mann and C. Holm, ESPResSo – an extensible simulation package for research on soft matter systems, Comput. Phys. Commun., 2006, 174(9), 704–727 CrossRef CAS.
- A. Arnold, O. Lenz, S. Kesselheim, R. Weeber, F. Fahrenberger, D. Roehm, P. Košovan and C. Holm, Espresso 3.1: Molecular dynamics software for coarse-grained models. in Meshfree Methods for Partial Differential Equations VI, ed. M. Griebel and M. A. Schweitzer, Lecture Notes in Computational Science and Engineering, Springer Berlin Heidelberg, 2013, vol. 89, pp. 1–23 Search PubMed.
- H. Tanaka, Y. Nishikawa and T. Koyama, Network-forming phase separation of colloidal suspensions, J. Phys.: Condens. Matter, 2005, 17, L143 CrossRef CAS.
- A. B. Dobroserdova and S. S. Kantorovich, Self-diffusion in monodisperse three-dimensional magnetic fluids by molecular dynamics simulations, J. Magn. Magn. Mater., 2017, 431, 176–179 CrossRef CAS.
- M. J. Stevens and G. S. Grest, Phase coexistence of a stockmayer fluid in an applied field, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 5976–5984 CrossRef CAS.
- J. E. Martin, E. Venturini, G. L. Gulley and J. Williamson, Using triaxial magnetic fields to create high susceptibility particle composites, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69(2), 021508 CrossRef PubMed.
- N. Osterman, I. Poberaj, J. Dobnikar, D. Frenkel, P. Ziherl and D. Babić, Field-induced self-assembly of suspended colloidal membranes, Phys. Rev. Lett., 2009, 103(22), 228301 CrossRef CAS PubMed.
- F. J. Maier and T. M. Fischer, Critical nucleation mesh-size of coarsening transient colloidal networks, Soft Matter, 2015, 12(612), 614–618 Search PubMed.
- T. Demirci, J. Teiser, T. Steinpilz, J. Landers, S. Salamon, H. Wende and G. Wurm, Is there a temperature limit in planet formation at 1000k?, Astrophys. J., 2017, 846(1), 48 CrossRef.

## Footnote |

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sm00796e |

This journal is © The Royal Society of Chemistry 2018 |