Time-dependent clustering and magnetization in magnetic colloidal suspensions

Abstract

This work analyzes the influence of the time-dependent clustering aggregation process on the transient and equilibrium magnetization of a monodisperse magnetic colloidal suspension under a uniform magnetic field via Brownian dynamics simulations. The clustering aggregation process is characterized by microstructural properties, such as the nucleation-growth factor, 〈n_c(t)〉, mean cluster size, 〈N_c(t)〉, kinetic exponent, z, effective radius, 〈R_eff〉, and radial distribution function, g(r). These are analyzed in terms of the volume fraction, ϕ, the dipolar coupling parameter, λ, and the Langevin parameter, α. Here, λ and α measure the magnetic dipole–dipole interaction energy and the magnetic field–dipole interaction energy relative to the thermal energy, respectively. The magnetization in transient and equilibrium regimes is analyzed relative to these microstructural properties for different values of ϕ, λ, and α. The analysis of the microstructural properties reveals a reduction in the dipolar chain growth at higher λ and ϕ values in the range of 1 < α < 10, which contrasts with the increase observed for low values of the same parameters. This reduction is caused by the lateral interactions between the chains formed. For higher ϕ and λ values, these interactions lead to side-by-side coupling of the long-dipolar chains that enhances the transient and equilibrium magnetization. The equilibrium magnetization values have been compared with some predictive models, showing a significant discrepancy at 0.01 ≤ α ≤ 10, which involves the aforementioned range of α. Also, the Langevin magnetic susceptibility, χ_L, used in these models provides a way to characterize dilute suspensions with strong magnetic interparticle interactions (χ_L ≥ 0.09). These results may contribute to formulating more accurate models to predict magnetization of these suspensions.

---

1 Introduction

A magnetic colloidal suspension is a dispersion of tiny magnetic particles immersed in a nonmagnetic carrier fluid.^1,2 Each particle comprises magnetic domains defined by the magnetic moments of individual atoms oriented in fixed directions. As the particle reduces its size on the nanometric scale (20–100 nm), these domains reduce to one characterized by an intrinsic permanent magnetic dipole moment,³m = M_sV_p, proportional to the volume of the magnetic core particle, V_p = πd³/6, where M_s is the saturation magnetization and d is the diameter of the particle.

Because of the nanometric length scale of the magnetic colloidal particles, they exhibit Brownian motion due to the thermal fluctuations of the carried fluid.^4–7 These prevent particles from settling out. A surfactant layer is applied on the particle surfaces to prevent agglomeration due to short-range attractive forces.^8–10

In a magnetic colloidal suspension, particles interact via their magnetic dipole moments, which induces the aggregation process of self-assembly, forming microstructures of different lengths and configurations.^11–15 This permits the tunability of these suspensions for use in different applications.^16–26 Furthermore, this cluster aggregation process and the cluster morphologies depend not only on the dipole–dipole interaction strengths but also other factors, such as an applied external magnetic field, particle concentration, bidispersity,²⁷etc.²⁸

Concerning the dipolar magnetic interaction, this is measured using the magnetic dipolar coupling constant, λ = μ₀m²/2πd³k_BT, which allows the clustering aggregation process when the anisotropic dipole–dipole interaction energy of two contacting particles, m²/d³, exceeds the thermal energy, k_BT, of the fluid, where k_B is the Boltzmann constant and T the absolute temperature. In the absence of an external magnetic field, rich cluster morphologies are formed according to the value λ.^15,29 As the dipolar interaction energy increases, there is a transition from separated magnetic particles³⁰ (λ < 1) to diverse structures. In this way, when the dipolar interaction energy is comparable to the thermal energy³¹ (λ = 1), there are a few pairs of particles (that is, dimers) in the suspension. However, when the dipolar interaction energy is much stronger than the thermal energy (λ ≫ 1),³² the latter cannot overcome the dipolar interactions. Therefore, the suspension develops different types of cluster morphologies, such as dimers, trimers, quadrumers, worm-like structures,³³ single flexible/semi-flexible chains,^{14,31,34–36} ring (flux-closure) structures,^11,12,37 and three- and four-fold junctions.³⁸

In the presence of an external magnetic field, the energy resulting from the interaction between the magnetic field and the magnetic dipoles, μ₀mH, competes with the thermal energy. The Langevin parameter, α = μ₀mH/k_BT, measures this competence, where H is the strength of the magnetic field. Depending on the field strength value (or the value α), the sizes and morphologies of the zero-field structures already formed (α = 0) are manipulated.^30,37,39 In addition, chain-like assembly structures are also induced by the field.^40–42 Several studies have shown that randomly distributed chains (at zero-field) form long dipolar chains oriented in the field direction for strong magnetic fields^31,37,41,43 (α ≫ 1) and strong magnetic dipole–dipole interactions (λ ≫ 1). Here, the short dipolar chains align in the field direction by forming ordered long dipolar chains or columnar structures.^41,43–45

In the case of particle concentration, it can be measured from the particle volume fraction, ϕ = nπd³/6, defined as the total volume of the N particles divided by the total volume of the colloidal suspension V, where n = N/V is the particle number density. In the limiting case of a very dilute suspension (ϕ ≪ 1), there is no evidence of any aggregates, even for strong magnetic dipolar interactions (λ ≫ 1). When the particle concentration increases, the places for the particles to move freely are reduced. They are then closed enough such that the translational Brownian diffusion is reduced and the magnetic dipole–dipole particle interaction is intensified. This effect on the dipolar interaction energy enhances the clustering aggregation process. Thus, for higher concentrations, it is known that short dipolar chains switch to longer chains to reach sufficient dimensions to be compared with the size of the colloidal system.^34,45,46

While the magnetic dipolar interaction (λ), the magnetic field (α), and the particle volume fraction (ϕ) affect the clustering aggregation process, they likewise influence macroscopic properties of the suspension, such as magnetization.^{33,42,45,47,48} This is one of the magnetic properties of fundamental and applicative interest in studying magnetic colloidal suspensions. Over the years, many theoretical, numerical, and experimental studies have been carried out to accurately predict magnetization.^{39,47,49–58} It has been found that, in the limit of low concentrations (ϕ ≪ 1), the magnetic suspensions behave as an ideal paramagnetic gas (negligible interparticle dipole–dipole interactions). Only the interaction of the external magnetic field with the individual magnetic dipoles plays a dominant role. Thus, following Langevin's classical theory, the equilibrium magnetization of these suspensions can be written as^2,59,60 〈M_eq〉 = ϕM_sL(α), where L(α) = coth(α) − 1/α is the Langevin function, and ϕM_s is the saturation magnetization of the suspension. However, when the magnetic interparticle dipole–dipole interactions become relevant, the Langevin magnetization is not a good fit to predict the equilibrium magnetization of the suspension correctly.^50–52,54 This occurs due to the formation of dipolar structures (e.g., dipolar chains)^{33,39,47,56,61,62} induced by the interparticle dipole–dipole interactions.

Consequently, several theoretical models for predicting magnetization have been developed, which consider in their framework the mechanisms of chain aggregation or dipole–dipole interactions.^{32,50–55,57,63–67} However, the discrepancy observed between the predicted magnetization and the data, especially when the interparticle magnetic interactions are strong, indicates that the influence of the clustering aggregation process on this property remains to be thoroughly investigated. This highlights the need for a rigorous specification and comprehensive description of the aggregation mechanisms to advance the models' accuracy. Therefore, the authors consider that an analysis of the temporal evolution of the clustering aggregation process may provide specific insights into the influence of this process on magnetization, which could contribute to improving predictive models.

The present work then analyzes the influence of the temporal evolution of the clustering aggregation process and cluster morphology on the transient and equilibrium magnetization of a magnetic colloidal suspension. By means of Brownian dynamics (BD) simulations, a qualitative and quantitative microstructural characterization is developed. The qualitative description is carried out through snapshots that show first-hand a direct observation of the microstructure at different instants of the initial nucleation and growth stages of the clustering aggregation process. Quantitatively, this latter is characterized by measuring microstructural properties such as the mean cluster size, 〈N_c(t)〉, the nucleation growth factor, 〈n_c(t)〉, the kinetic exponent, z, the effective radius, 〈R_eff〉, and the radial distribution function g(r) at a distance r from a reference particle for different values of magnetic field strength (or α), magnetic dipolar strength (or λ), and particle volume fraction (ϕ). In addition, to independently assess the influence of these parameters on the clustering aggregation process and magnetization, simulations were performed in the dilute regime.

The sections in this work are organized as follows: the model of the magnetic colloidal suspension is described in Section 2. The governing equations that describe the dynamics of the translational and rotational motion of the magnetic colloidal particles are presented under the BD methodology in Section 3. The definitions of the microstructural properties and magnetization are provided in Section 4. The microstructure developed for set values of λ, α, and ϕ, the results of the microstructural properties, the time-dependent magnetization, and the equilibrium magnetization are presented in Section 5. Finally, the main conclusions of this work are presented in Section 6.

2 Model system

The model system consists of a monodisperse suspension of N identical rigid spherical magnetic colloidal particles of radius a immersed in an incompressible Newtonian fluid of dynamic viscosity η and density ρ. This carrier fluid is quiescent in the laboratory reference frame x–y–z, as shown in Fig. 1. Each particle bears a permanent magnetic dipole moment m = mu of magnitude m and orientation u, rigidly fixed at its center of mass, oriented along the x′-axis of the rigid body reference frame x′–y′–z′ on the particle. This model also considers a uniform external magnetic field, H, applied along the z-axis of the laboratory reference frame x–y–z, as illustrated in Fig. 1. Each particle also has a thin steric layer of thickness δ to prevent particle overlapping in suspension dynamics.

Fig. 1 Monodisperse suspension of magnetic interacting colloidal particles of radius a under a uniform magnetic field H. The permanent point magnetic dipole moment m of each particle is fixed at its center of mass. Each dipole moment direction is established by a black arrow into the particle. The dark red and light blue part of each particle represents its north and south pole, respectively. A steric layer of thickness δ is covering each particle.

3 Brownian dynamics simulations

The starting point of the BD simulation method is the Langevin equations^68–72 that model the force and torque balance that act on the colloidal particles. In these equations, the terms related to the inertia of the particles are negligible because the numerical resolution timescale to solve the equations is significantly larger than the translational and rotational particle inertial relaxation times. Thus, for a colloidal suspension composed of rigid interacting magnetic particles, in the presence of a uniform magnetic field, the Langevin equations for translation and rotation read as follows:⁷³where the subscripts i and j denote the generic particles of the suspension. Here, F^H_i and T^H_i are the hydrodynamic viscous force and torque acting on the particle i, respectively; F^M_i and T^M_i are the force and torque exerted by the external uniform magnetic field on the particle i, respectively; F^D–D_ij and T^D–D_ij are the magnetic dipole–dipole interaction force and torque, respectively, exerted on the particle i by a particle j. F^R_ij is the repulsive force exerted by a particle j on the particle i when their steric layers are overlapped. F^B_i and T^B_i are the Brownian force and torque, respectively, acting on the particle i.

Due to the very small size of the particles, the flow surrounding each one is the well-known Stokes flow. This means that the hydrodynamic Stokes drag force, F^H_i, and torque, T^H_i, for a spherical particle i, with hydrodynamic radius a, moving with linear, U_i, and angular, Ω_i, velocities in a Newtonian fluid, are calculated as follows:^74,75

F^H_i = −R_FU·U_i, (3)

T^H_i = −R_TΩ·Ω_i, (4)

where R_FU = 6πηaI and R_TΩ = 8πηa³I are the hydrodynamic resistance tensors, for translation and rotation, in the absence of hydrodynamic interactions; I is the unit isotropic tensor. Here, many-body and lubrication hydrodynamic interactions are neglected due to the dilute regime assumed to be the basis of the suspensions analyzed in this work.

The magnetic dipole–dipole interaction force and torque exerted by the particle j on the particle i are obtained from the dipolar interaction potential:⁷⁶

where m_i and m_j are the magnetic dipole moments of the particles i and j, as shown in Fig. 2, respectively; r_ij = r_i − r_j is the separation vector between the positions of the particles i and j, r_ij = |r_ij| is the magnitude of the vector r_ij (center-to-center distance between the particles) and μ₀ = 4π × 10⁻⁷ N A⁻² is the magnetic permeability of free space.

Fig. 2 Steric repulsion force, F^R_ij, on particle i exerted by particle j when their steric layers, with thickness δ, are overlapped. The radius of each particle is a. The vector r_ij is the center-to-center vector from particle j to particle i. Both particles i and j exhibit centered point dipole magnetic moments m_i and m_j, respectively.

The dipole–dipole magnetic interaction force and torque between the particles are then obtained by applying the translational, ∇ = ∂/∂r, and orientational, gradient operators, respectively, to the magnetic potential (5), as follows:

F^D–D_ij = −∇Φ^D–D_ij(r_ij, m_i, m_j), (6)

Thus, the explicit equations for the magnetic dipole–dipole interaction force and torque, respectively, read⁷³

where r̂_ij = r_ij/r_ij is the unit vector of r_ij.

The resulting magnetic dipolar interaction force (8) and torque (9) are coupled. The magnetic dipole–dipole torque induces a rotation of the particle point-dipole moments. This causes a change in the orientations of the dipoles, which in turn affects the dipolar force between the particles, according to (8). Similarly, this interparticle dipole–dipole force induces a translation of the particle point-dipole moments. The relative position between the particles then affects the relative orientation of the point-dipole moments, which affects the dipolar torque between the particles, according to (9), and so on.

The field–particle interaction between the uniform magnetic field, H, and the point–dipole moment, m_i, induces a magnetic torque on the particles, given by

where Φ^M(m_i, H) = −μ₀m_i·H is the external magnetic potential on the ith-particle. This external potential also induces a magnetic force on the particles. But in this case it is zero, F^M_i = 0, because the external magnetic field is uniform and m_i is constant.

The steric repulsion energy between the magnetic particles i and j, as illustrated in Fig. 2 is¹⁰

where λ_R is called the steric layer factor, δ is the thickness of the steric layer, and T is the absolute temperature of the suspension.

Then, the repulsive force to prevent particles from approaching so close is calculated using F^R_ij = −∇Φ^R_ij(r_ij), which in turn reads as

Brownian motion is due to the thermal fluctuations of fluid molecules acting on suspended particles through a stochastic force and torque, called Brownian force and torque, characterized by the fluctuation–dissipation theorem^77–79 as follows

〈F^B_i〉 = 0 and 〈F^B_i(0)F^B_i(t)〉 = 2k_BTR_FUδ(t), (13)

〈T^B_i〉 = 0 and 〈T^B_i(0)T^B_i(t)〉 = 2k_BTR_TWδ(t), (14)

where δ(t) is the Dirac delta function.

The translational and rotational evolution equations of the ith particle are obtained from (1) and (2), and read^80,81

〈Δr̃^B_i〉 = 0 and 〈Δr̃^B_iΔr̃^B_i〉 = 2IΔt̃, (15b)

where each magnitude with a tilde (^∼) denotes a dimensionless quantity. Here, Δr̃_i and Δ[small theta, Greek, tilde]_i are the changes in the linear and angular positions of the particles during the time step Δt̃; Δr̃^B_i(Δt̃) and Δ[small theta, Greek, tilde]^B_i(Δt̃) are the random linear and angular displacements due to Brownian motion that have zero means and covariances, 2D^traI for translation and 2D^rotI for rotation, respectively, where D^tra = k_BT/6πηa and D^rot = k_BT/8πηa³ are the translational and rotational diffusion coefficients of a single isolated particle. In (15a) and (15c), r̃_i has been nondimensionalized by the characteristic particle size a, and the time by the translational diffusion time scale, τ_B = a²/D^tra.

In (15a) and (15c), λ and α are dimensionless parameters that characterize the mechanisms of formation and manipulation of the suspension microstructure. These parameters have already been introduced in Section 1; bear in mind that λ is the dipolar coupling constant parameter that measures the magnetic interparticle dipole–dipole interaction energy relative to the thermal energy, and α is the Langevin parameter that measures the field–particle interaction energy relative to the thermal energy.

In addition, λ and α can be interpreted as dimensionless forms of the magnetic dipolar and magnetic field strengths, respectively. Another relevant physical property considered in this work is the particle volume fraction, ϕ, also introduced in Section 1, which determines the suspension concentration. Recalling and rewriting these parameters here can be expressed as follows:

where ñ stands for the dimensionless form of the particle number density n.

To establish a rigid body reference frame x′–y′–z′ where the particle magnetic moments have a permanent magnetization (i.e., fixed orientation) along the x′-axis as illustrated in Fig. 1, the laboratory reference frame x–y–z suffers a transformation (or three successive angular rotations) characterized by the well-known Euler angles θ, ϕ, and ψ.

These rotations are measured via a transformation matrix A in terms of Euler angles. Adopting an x-convention (that is, second rotation about the intermediate x axis) in the rotations and expressing the Euler angles in terms of the quaternions e₀, e₁, e₂, and e₃ (also known as Euler parameters) in order to free numerical computations from singularities, the transformation matrix reads⁸²

where the quaternion parameters satisfy the relation e₀² + e₁² + e₂² + e₃² = 1.

The transformation matrix A is then used to transform the magnetic moment of each particle from the rigid-body reference, to the laboratory reference frame, as follows:

where A^T is the transpose matrix of A and is permanently oriented along the x′-axis with the unit vector i_x′.

The time evolution of the quaternion parameters is calculated by first evaluating in an instant of time the change in the angular displacement of the particle magnetic moments in the rigid laboratory reference frame, as follows:

where is the dimensionless magnetic torque; note that according to our model system, u′ = i_x′ and H̃′ = A·H̃ with H̃ = i_z along the z-axis, as illustrated in Fig. 1.

The angular displacement obtained from (19) is then used in the matrix eqn (20) to determine the change of the quaternion parameters at each time step:^83,84

where the quaternion of each particle is normalized in each time step:e = (e₀, e₁, e₂, e₃) and Δe = (Δe₀, Δe₁, Δe₂, Δe₃).

4 Microstructural and magnetic properties

4.1 Microstructural properties

In order to achieve a good characterization of the temporal behavior of the microstructure suspensions, several microstructural properties are measured over time. These properties are the nucleation-growth factor 〈n_c(t̃)〉, the mean cluster size 〈N_c(t̃)〉, the kinetic exponent, z, the effective radius 〈R_eff〉, and the radial distribution function, g(r).

In this study, two particles i and j aggregate with each other when r_ij/2a < 1.2.⁸⁵ The particles that do not form clusters are called singlets and are also considered as clusters of the suspension.

The nucleation-growth factor, 〈n_c(t̃)〉, is defined as the ratio of the number of clusters, regardless of the number of singlets N_singlets, and the sum of the number of particles per cluster in the suspension. Regarding k three-dimensional clusters formed in the colloidal suspension (see Fig. 3), 〈n_c(t̃)〉 is defined as

where the angular brackets 〈〉 represent an ensemble average on the total of the clusters, and N_ci stands for the number of particles making up the ith cluster (c_i) of the suspension.

Fig. 3 Outlining of k three-dimensional clusters formed by magnetic particles of a monodisperse magnetic colloidal suspension due to interparticle dipolar interactions. Each cluster c_i comprises N_ci particles and is enclosed in the smallest possible sphere with radius R_ci, where i = 1, 2, …, k.

The mean cluster size, 〈N_c(t̃)〉, indicates the mean number of particles per cluster formed as a function of time, considering singlets. This microstructural property is calculated as follows

Note that the mean cluster size at long times asymptotically behaves like a power-law function^86–88 as a function of time as

〈N_c(t̃)〉 ∼ t̃^z, (24)

where z is called the kinetic exponent. In addition, this work examines the influence of the magnetic field on the mean cluster size and the z exponent values.

The effective radius 〈R_eff〉 of a suspension is the weighted average of the radii of the smallest possible spheres that surround the formed clusters, as shown in Fig. 3:

where R_ci is the radius of the smallest sphere that encloses the number of particles per cluster, N_ci.

Finally, the radial distribution function, g(r), depicts the probability of finding a particle at a distance r from a reference particle, compared to the same probability expected for a completely random distribution.

In dimensionless form, this property is calculated as follows:

where Δr̃ is the width of the dimensionless bin and r̃_ij the dimensionless center-to-center distance between particles i and j.

4.2 Macroscopic property: magnetization

Instantaneous and equilibrium magnetization are the macroscopic magnetic properties analyzed in this work. The instantaneous (or time-dependent) magnetization is obtained by averaging the particle magnetic moments of the suspension via⁸⁹where m̃_i(t̃) can be obtained from (18). Here, the angle brackets denote an ensemble average on the number of particles, and the magnetization is dimensionless via the saturation magnetization of the suspension, ϕM_s.

For magnetic dipole moments permanently magnetized along the x′-axis where m̃′(t̃) = i_x′ + 0i_y′ + 0i_z′, the instantaneous magnetization parallel to the field direction (or z-axis) reads

which is the only magnetization reported in this work, hereafter denoted by 〈M̃(t̃)〉 = 〈M_z(t̃)〉/ϕM_s. The other components of the magnetization perpendicular to the magnetic field,become zero for each time step, that is, 〈M̃_x(t̃)〉 = 〈M̃_y(t̃)〉 = 0.

The equilibrium magnetization of the suspension, 〈M̃_eq〉 = 〈M_eq〉/ϕM_s, is then obtained by evaluating (28) at times long enough such that the instantaneous magnetization achieves its saturation (or steady-state) regime. It can be expressed as

In this work, the equilibrium values of the magnetization are obtained by averaging the last 10% of the magnetization data. This corresponds to the saturation tendency of the data on linear scale in this range. The equilibrium magnetization values obtained for different values of field strength (α), dipolar coupling strength (λ), and volume fractions (ϕ) are then compared with Langevin magnetization.

5 Results

The number of particles used in the simulations is N = 500. For particles N ≥ 500, our system behaves statistically independently. The initial configuration is such that the initial positions and orientations of the particles are randomly distributed in a cubic simulation box. The simulations are carried out in the simulation box with periodic boundary conditions where the minimum image conventions are adapted.⁶⁸ The BD simulations implement the steric repulsion potential to prevent particle overlapping and agglomeration. Here, the dimensionless thickness of the steric layer that covers each particle is [small delta, Greek, tilde] = δ/a = 0.01 with the steric layer factor λ_R = 1300.⁹⁰

Three different values of ϕ have been selected in order to analyze the influence of this parameter on microstructure and magnetization. In addition, the parameter λ ranges from 5 to 75 in order to analyze the effect of weak and strong magnetic dipole–dipole interactions between the particles. Finally, the parameter α ranges from 0 to 1000 to characterize the application of a magnetic field from weak to strong intensities. Table 1 shows the values of the dimensionless parameters of these parameters used to perform the simulations. For each set of ϕ–λ–α values, ten repetitions are performed to establish good statistics in the results. The end time of the simulations is 1000 with a time step Δt̃ = 10⁻⁴ dimensionless units.

Table 1 Set of values of particle volume fraction (ϕ), dipolar coupling parameter (λ), and Langevin parameter (α) for simulations of 500 particles

ϕ	λ	α
10^−4	5	0.01
	15	0.03
	45	0.1
	75	0.25
		0.5
5 × 10^−4	5	1
	15	1.5
	45	2.5
	75	5
		10
10^−3	5	25
	15	50
	45	100
	75	400
		1000

---

5.1 Microstructural characterization of the clustering process

5.1.1 Initial nucleation and growth stages. The temporal evolution of the clustering aggregation process passes through two well-defined characteristic stages: the initial nucleation stage and growth (or coalescence) stage.^85,91,92 The first stage initiates the formation of several short clusters, evolving from singlets (monomers) to short-dipolar chains such as doublets (dimers), triplets (trimers) and quadruplets (quadrumers), as shown in Fig. 4. As time goes on, the growth stage is initiated. This consists of the aggregation of the short clusters, leading to larger structures, frequently long chain-like structures, that continue growing, as also illustrated in Fig. 4.

Fig. 4 Nucleation and growth stages of the clustering aggregation process of a monodisperse magnetic colloidal suspension for ϕ = 10⁻³, λ = 15, and α = 10. The vertical black arrows represent the magnetic field, H, applied along the z-axis.

The formation of these dipolar chains is influenced by the interplay between Brownian motion and the magnetic dipole–dipole interaction, measured by λ. Hence, for weak dipolar strengths, such as λ = 5, thermal fluctuations inhibit the clustering aggregation process (see Fig. 12) due to the predominant Brownian diffusivity of singlets. In contrast, when dipolar strength predominates in Brownian motion (here, λ = 15, λ = 45, and λ = 75), the imminent formation of clusters in the nucleation and growth stages is enhanced as λ increases, as shown in Fig. 4, 8 and 13.

The particle volume fraction, ϕ, also affects the clustering aggregation process. The population of short-dipolar clusters in the nucleation stage and the size of the dipolar chains in the growth stage increase as ϕ increases. The presence of an external magnetic field influences the ordering of the dipolar chains in the nucleation and growth stages. As the field strength increases, the long dipolar chains transform from randomly disordered and flexible to oriented rigid dipolar chains along the field direction, as illustrated in Fig. 8.

In addition to showing a direct observation of the temporal evolution of the clustering aggregation process, a quantitative description is also given in terms of the nucleation-growth factor, 〈n_c(t̃〉), the mean cluster size 〈N_c(t̃)〉, the kinetic exponent z, the effective radius, 〈R_eff〉, and the radial distribution function, g(r̃), as illustrated in Fig. 5, 6, 7, 9 and 10, respectively.

5.1.2 Nucleation-growth factor. The nucleation and growth stages are related to the number of short-dipolar and long-dipolar (chain) clusters formed. The nucleation stage leads to an increase in the population of clusters over time, which decreases in the growth stage because clusters aggregate to form larger structures. The use of the nucleation-growth factor, 〈n_c(t̃)〉, allows a quantitative description of both stages and the influences of λ, ϕ and α on them.

The temporal evolution of 〈n_c(t̃)〉, illustrated by the curves in Fig. 5, shows the initial increasing behavior of 〈n_c(t̃)〉 due to the increase of the number of clusters, corresponding to the nucleation stage. Afterward, the decay of the 〈n_c(t̃)〉 curves characterizes the aggregation of these clusters, causing a decrease in their population, which indicates the growth stage.

Fig. 5 Temporal evolution of the nucleation-growth factor, 〈n_c(t̃)〉, for moderate (λ = 15) and very strong (λ = 75) magnetic dipolar strengths under very weak (α = 0.01), moderate (α = 1), and very strong (α = 1000) magnetic field strength values. In the figure, (a) and (c) correspond to λ = 15, whereas (b) and (d) correspond to λ = 75. The particle volume fraction varies from ϕ = 10⁻⁴ (a and b) to ϕ = 10⁻³ (c and d). N.S. stands for nucleation stage and G.S. for growth stage.

Since the minimum possible cluster formed with two or more particles is a doublet, the theoretical threshold of 〈n_c(t̃)〉 according to (22) is 0.5. This implies that the entire cluster population comprises only doublets. However, since 〈n_c(t̃)〉 never reaches this value, the cluster population in the growth stage comprises a variety of clusters that involve singlets, doublets, triplets, etc.

As a starting point, the effect of λ on 〈n_c(t̃)〉 is analyzed. The absence of a clustering aggregation process, as occurs for λ = 5 (see Fig. 12), results in 〈n_c(t̃)〉 being a constant zero. Then, the increase of λ promotes the clustering aggregation process. For low values, such as λ = 15, Brownian motion prevails and prevents particles from easily aggregating. Conversely, high λ values, such as λ = 75, ensure that random moving particles are easier to couple via a head-to-tail configuration. Thus, increasing λ (from λ = 15 to λ = 75) facilitates the aggregation process, which induces a greater number of short-dipolar clusters (〈n_c(t̃)〉) in the nucleation stage and a lower number of them (〈n_c(t̃)〉) in the growth stage, as observed in Fig. 5(a)–(d). In addition, the duration of the nucleation stage is also reduced, revealing that the clustering process develops faster with increasing λ.

In the same way, ϕ also affects the clustering aggregation process. Low ϕ values cause the particles to be far away from each other. This reduces the number of dipole–dipole interactions between particles, resulting in slow aggregation. In contrast, high values lead particles to be closer to each other, which enhances dipole–dipole interactions. Thus, the effect of ϕ on 〈n_c(t̃)〉 resembles the previous case, as observed in Fig. 5(a)–(c) and (b)–(d). In any case, the duration of the nucleation stage is reduced as ϕ increases from ϕ = 10⁻⁴ to ϕ = 10⁻³.

The influence of α on 〈n_c(t̃)〉 is also analyzed. Values of α ≠ 0 imply the presence of a magnetic field on the suspension. Singlets and clusters interact with the field due to the particle magnetic dipole moments, which tend to align along the field direction according to the magnitude. However, this trend competes against thermal fluctuations: the stronger/weaker the magnetic field, the stronger/weaker the trend compared to thermal fluctuations.⁹³

For low α values, thermal fluctuations have a more significant effect on the particle orientations than the magnetic field. The random rotation of particles makes it difficult to achieve a head-to-tail configuration, and consequently, the formation of clusters. In contrast, for high α values, the magnetic field overcomes the thermal fluctuations. However, the particle dipole moments still exhibit random rotation, but are restricted to rotation around the magnetic field direction. In this case, dipole moments tend to align along the field direction, which favors particles achieving the head-to-tail configuration.

Very high values of α make all magnetic dipole orientations parallel, hardly enhancing the aggregation process. Faster aggregation is expected as α increases, which implies greater 〈n_c(t̃)〉 in the nucleation stage and lower 〈n_c(t̃)〉 in the growth stages.

However, expected behavior occurs only at α > 1 for very low ϕ and λ values, as observed in Fig. 5(a). Here, particles are very far from each other and are weakly interacting. Therefore, similar orientations due to their high α drive aggregation, and hence the formation of more clusters. For α < 1, the behavior of 〈n_c(t̃)〉 differs from what is expected. Here, 〈n_c(t̃)〉 is independent of the field effect, and one curve (α = 1) overlaps the other (α = 0.01), as observed in Fig. 5.

Furthermore, as ϕ and λ increase, the gradual collapse of all 〈n_c(t̃)〉 curves into one, as observed in Fig. 5(b)–(d), suggests independence of the clustering aggregation process from α. This seems to be supported by the unchangeable duration of the nucleation stage in each system, whereas α takes greater values.

Nonetheless, the 〈n_c(t̃)〉 curves for α > 1 exhibit a slight shift to the right as well as ϕ and λ reaching very high values (see Fig. 5(d)). This leads to a lower number of clusters in the nucleation stage and a higher number in the growth stage, compared to those in the low α regime. This suggests a possible delay of the clustering aggregation process for α > 1, even when the magnetic dipole moments are aligned with the field. Therefore, α still has an effect on the clustering process despite all the 〈n_c(t̃)〉 curves seeming to collapse into one.

The effect of α on the clustering aggregation process should be analyzed. In order to gain insight into this and complete the description of the clustering process, it is necessary to analyze the cluster size and morphology via the mean cluster size and the effective radius.

5.1.3 Mean cluster size and kinetic exponent. Due to the progressive aggregation between particles and clusters, the cluster sizes constantly change over time. To get an overview of these temporal changes, the transient behavior of the mean cluster size, 〈N_c(t̃)〉, is analyzed. The absence of a plateau behavior of the resulting 〈N_c(t̃)〉 curves shown in Fig. 6 suggests that the clustering aggregation process continues to evolve after t̃ = 1000.

Fig. 6 Influence of the magnetic dipolar strength (λ = 15, 45, 75), particle volume fraction (ϕ = 10⁻⁴, 5 × 10⁻⁴, 10⁻³), and magnetic field strength (α = 0.01, 1, 1000) on the mean cluster size 〈N_c(t̃)〉 as a function of time (t̃).

In order to understand the influences of λ and ϕ on 〈N_c(t̃)〉, it is necessary to analyze the force and torque of the dipole–dipole magnetic interaction. Both drive the clustering aggregation process, inducing particles to achieve the head-to-tail configuration. The magnetic torque rotates the interacting dipoles to reach a parallel orientation, whereas the force sticks them to each other in the head-to-tail configuration. Consequently, increasing both the force and torque, this configuration is achieved faster and by more particles, which in turn enhances the clustering process. Then, according to (8) and (9), there are two ways to increase this force and torque.

The first way is by increasing the magnitude of the dipole moments, which implies an increase of the corresponding λ. This strengthens the magnetic dipole–dipole interactions over thermal fluctuations. Then, more particles tend to aggregate in the same cluster, which implies an increase in cluster size. Therefore, increasing λ causes higher 〈N_c(t̃)〉 values over time, as observed in Fig. 6.

The second way is by having the particles closer to each other. This is possible by increasing ϕ. The force and torque between two particles decay with r⁴ and r³, respectively. Thus, if the particles are closer to each other (higher ϕ), both force and torque hardly increase. Thus, the increase in ϕ also leads to the increase in 〈N_c(t̃)〉, as observed in Fig. 6.

In addition, the 〈N_c(t̃)〉 curves exhibit distinctive behavior when α varies. The enhancement triggered by the increase of α in the alignment of the magnetic dipole moment of the particle must contribute to faster aggregation and, therefore, to higher cluster sizes. Thus, it is expected that the 〈N_c(t̃)〉 curves arrange themselves one on top of the other as α increases.

The collapse of the 〈N_c(t̃)〉 curves in the low α regime indicates the independence of cluster sizes from the magnetic field effect. This confirms the predominance of thermal fluctuations over the magnetic field in this regime. For α > 1, expected behavior is observed, but only for λ = 15, disappearing as ϕ increases, as shown in Fig. 6. When increasing λ, the collapse of the 〈N_c(t̃)〉 curves at ϕ = 10⁻⁴ agrees with the behavior of the corresponding 〈n_c(t̃)〉 curves. However, as λ and ϕ reach higher values, the 〈N_c(t̃)〉 curves for α > 1 exhibit a drop below those for the low α regime. This comes from the slight shift of the 〈n_c(t̃)〉 curves shown in Fig. 5(d): the lower/greater number of clusters in the nucleation/growth stage implies smaller clusters in the system. Therefore, the 〈N_c(t̃)〉 curves decay at α > 1.

In order to better justify the description given above, the behavior of 〈N_c(t̃)〉 for several values of α is analyzed. Starting from a certain time, all the 〈N_c(t̃)〉 curves display exponential growth. Applying curve fitting to this exponential behavior, the kinetic exponent z of each curve is obtained. This is a characteristic parameter of the aggregation rate, the value of which depends on different characteristics such as the type of magnetic particle, the applied magnetic field, etc.^{32,86,87,94–96}

Then, z represents the temporal behavior of 〈N_c(t̃)〉. In support of this, z exhibits an increasing behavior, as observed in Fig. 7, with higher values of ϕ and λ, reflecting the corresponding growth of 〈N_c(t̃)〉 under the same conditions, as shown in Fig. 6. This is expected compared to the behaviors of z observed in similar systems studied.^65,97 Furthermore, z shows a fairly constant behaviour in the low α regime (see Fig. 7), which agrees with the collapse of all the 〈N_c(t̃)〉 curves in this regime, as observed in Fig. 6.

Fig. 7 Influence of uniform magnetic field (α values) on the kinetic exponent, z, of a magnetic colloidal suspension. Keeping ϕ constant, z is compared at λ = 15, 45 and 75 under different values of α.

However, when the value of α varies while ϕ and λ remain constant, z shows different behaviors. In this way, if λ = 15, z increases as α > 1, which corresponds to the increase of the 〈N_c(t̃)〉 curves in the same range of α. This reveals that the application of a magnetic field under the aforementioned conditions contributes to the growth of the clusters, regardless of the value of ϕ. This behavior differs from that observed in paramagnetic particles.⁸⁸ This is due to the constant and stronger magnetic dipole moment of the particles, which induces the clustering aggregation process despite the corresponding low λ value.

In contrast, for higher λ values, such as λ = 45 and λ = 75, the exponent z experiences a constant behavior in very dilute suspensions (e.g. ϕ = 10⁻⁴) or undergoes a decrease from the low α regime onward, to reach a near-constant value at higher α, as observed in Fig. 7. The decay of z has also been observed in previous works, for superparamagnetic particles and concentrated suspensions,^65,94 with values significantly lower than those reported in this work. This drop in z also confirms the decay of the 〈N_c(t̃)〉 curves in Fig. 6 for α > 1, mainly in the range 1 < α < 10. This interval will be relevant later. Beyond this, the 〈N_c(t̃)〉 curves collapse onto a single one due to the constant behavior of z. This is because all magnetic dipole moments are kept parallel to each other and, therefore, the increase in α does not perform an enhancement on the clustering aggregation process.

In general, the kinetic exponent z reveals differences in the clustering aggregation process for the suspensions studied, with the magnetic field inducing an enhancement of this process in some cases and slowing the process for other ones.

Hence, the effect of ϕ and λ on the clustering aggregation process is to reach larger dipolar clusters through a faster aggregation as these parameters increase. In contrast, the effect of α is similar, but only for very low values of ϕ (ϕ = 10⁻⁴) and λ (λ = 15). For high values of these parameters, the increase of α beyond the low α regime slows down the clustering process. It is necessary to develop a cluster morphology analysis in order to understand this clustering process delay.

5.1.4 Microstructure morphology and effective radius. The most common cluster in the systems is the chain-like structure, as shown in Fig. 8. However, this can undergo distortions due to thermal fluctuations.⁹³ Because α stands for the competence between magnetic field alignment trend and thermal fluctuations, α must hardly influence the morphology of clusters.³⁷

Fig. 8 Microstructure in a magnetic colloidal suspension with a concentration of ϕ = 10⁻³ at t̃ ≈ 1000. From left to right, the dipolar coupling parameter, λ, increases (λ = 15, 45, 75). From bottom to top, the uniform magnetic field H, characterized by α, takes different values: α = 0, 0.01, 1 and 1000. Each arrangement of black arrows stands for the direction of the applied magnetic field.

These distortions produce a bending effect on chains when thermal fluctuations predominate (low α regime). This can enclose them, producing ring-like structures¹² (see Fig. 8). However, in the case of the predominance of thermal fluctuations over magnetic dipole–dipole interactions, as occurs for λ = 15, the rings adopt short-lived behavior because of their constant breakage caused by thermal fluctuations. Meanwhile, high values of λ lead to longer-lived rings due to the strong bonds produced by the magnetic dipole–dipole interactions, as occurs for λ = 45 and λ = 75.

As α approaches 1, the growing magnetic-field alignment trend leads to a gradual restriction of the bending effect, which prevents the formation of ring-like structures. Consequently, the flexible chains tend to spread out along the magnetic field direction, as shown for α = 1 in Fig. 8. For α > 1, all chains align with the magnetic field, whereas thermal fluctuations still produce misalignment along them. For very high α values, such as α = 1000, these misalignments become very slight on the chains, which adopt a rigid chain-like behavior.

To get an overview of the microstructure morphology, the effective radius, 〈R_eff〉, is calculated as a function of 〈N_c(t̃)〉 (25), and compared to the rigid chain-like behavior (〈R_effchains〉 = 〈N_c(t̃)〉). This comparison reveals two key points. The first concerns the initial behavior of the 〈R_eff〉 curves. According to Fig. 9 and the results discussed previously, the clusters behave on average like rigid chains throughout the nucleation stage. Then, the clustering aggregation process at this stage can be modeled as aggregation of short rigid chains. This important feature will be useful for the analysis of the transient behavior of magnetization.

Fig. 9 Relationship between effective radius, 〈R_eff〉, and mean cluster size, 〈N_c〉, for magnetic colloidal suspensions under different uniform magnetic fields (α = 0.01, 1, 1000). The systems have a volume fraction of ϕ = 10⁻³ and (a) λ = 15, (b) λ = 75. The solid line shows the theoretical 〈R_effchains〉 relationship between the same properties for a rigid chain. The snapshots show the common morphologies of the clusters in each system.

The second key point relates to the microstructure morphology at long times. In the low α regime, such as α = 0.01, the divergence of the 〈R_eff〉 curves from the rigid chain-like behavior (see Fig. 9) relies on the predominance of thermal fluctuations, leading to flexible chains and the formation of ring-like structures. However, this divergence occurs later for λ = 75 (Fig. 9(b)) than for λ = 15 (Fig. 9(a)). This suggests that increasing λ produces more rigid particle–particle bonds, which allows short clusters to conserve the rigid chain-like behavior for larger values 〈N_c〉. On the other hand, for very high values of α, such as α = 1000, the clusters behave like rigid chains. However, since ϕ and λ have high values (Fig. 9(b)), clusters tend to move slightly away from rigid chain-like behavior as they reach large sizes (〈N_c〉). This is a consequence of the presence of node-like bonds along some chains, as shown in the corresponding snapshot in Fig. 9(b).

In a different vein, as also observed in Fig. 9(b), the case of α ≈ 1 with ϕ = 10⁻³ and high λ induces the formation of a double chain-like structure. This is also evidenced in Fig. 10. Although the increase in g(r̃) at r̃ ≈ 2.0 as λ increases from 15 to 75 is mainly due to the decay of the singlet population, the formation of double chains also contributes to this. In the same way, these structures also account for the prominent increase in the peaks at values of r̃ around 4.0 and 6.0, as these indicate a significant increase in the population of particles at such distances, which is directly associated with the formation of double chains.

Fig. 10 Radial distribution function, g(r̃), of a monodisperse magnetic colloidal suspension for ϕ = 10⁻³, α = 1 and λ = 15 (solid red curve) or λ = 75 (dashed blue curve) at t̃ = 1000.

The formation of these clusters is related to the flexibility of the chain and its very long size. As shown in Fig. 11, for the low α regime, such as 0.01, this flexibility allows long chains to aggregate by sticking their ends together (ESI,† Video 1). This mechanism also works for α = 1, but the flexibility is restricted due to the magnetic field alignment trend. This hinders their ends from encountering rapidly and sticking with each other. However, a new mechanism is developed that is possible if long chains exhibit flexibility (ESI,† Video 2).⁹⁸ As shown in Fig. 11, the chains are bent, laterally sticking and closing like a zipper. Otherwise, the lack of flexibility makes two long rigid chains interact laterally but not form a double chain, as observed for α = 1000 in Fig. 11 (ESI,† Video 3).

Fig. 11 Cluster–cluster magnetic interaction in a monodisperse suspension with ϕ = 10⁻³, λ = 75 and varying α from 0.01 to 1000.

Consequently, Fig. 11 explains the decay of the 〈N_c(t̃)〉 curves and the slight shift to the right of the 〈n_c(t̃)〉 curves at high ϕ and λ. For the low α regime, such as 0.01 and 1, the clusters aggregate with each other by joining their ends or forming double chains. This causes similar aggregation rates, which leads to overlapping of the corresponding 〈n_c(t̃)〉 and 〈N_c(t̃)〉 curves. However, at high α values, such as α = 1000, the unsuccessful lateral interaction between long rigid chains delays chain aggregation, which leads to a lower/greater number of clusters in the nucleation/growth stage (slight shift to the right) and a delay in the growth of mean cluster size (decay of the 〈N_c(t̃)〉 curve).

In general, the effect of α on the clustering process is weakly related to the increase in the number of clusters and cluster sizes, but is hardly related to the microstructure morphology, which in turn affect the aggregation rate of a suspension. This completes the description of the influence exerted by ϕ, λ, and α on the temporal evolution of the suspension microstructure. Using this information, the behavior of the magnetization of the suspensions will be explained.

5.2 Magnetization

5.2.1 Clustering aggregation process effect on the magnetization. The Langevin function shows the influence of α on 〈M̃_eq〉 of a suspension without the clustering process. According to this, the stronger the magnetic field (higher α), the greater the 〈M̃_eq〉 to reach the saturation magnetization value, ϕM_s. This occurs because of the progressive improved alignment of the particle dipolar moments with the magnetic field direction as α increases, to achieve perfect alignment.

According to previous works,⁹⁹ the clustering process does not occur from λ = 1. This agrees with the simulations made in this study for λ = 5 (see the frames in Fig. 12). Then, due to the absence of a clustering process, the Langevin function must accurately predict 〈M̃_eq〉. This can be verified by tracking the transient behavior of 〈M̃(t̃)〉 (29) to reach the corresponding 〈M̃_eq〉 value, which results from the long-term balance between all the effects acting on particle magnetic dipole moment orientations. As shown in Fig. 12, 〈M̃_eq〉 is achieved for any α value. The comparison between the obtained 〈M̃_eq〉 values with those predicted by Langevin function will be made later.

Fig. 12 Time evolution of the magnetization, 〈M̃(t̃)〉, of a magnetic colloidal suspension for ϕ = 10⁻³ and λ = 5 under different magnetic field strengths (α values). Snapshots at t̃ = 1000 reveal the absence of the clustering aggregation process in these suspensions.

The absence of the clustering aggregation process for λ = 5 implies that the long-term balance is performed only between thermal fluctuations and magnetic field–particle interactions. Then, following Fig. 12, this balance is achieved in a very short time interval, since 〈M̃(t̃)〉 quickly reaches its 〈M̃_eq〉 value.

In contrast, the development of the clustering process results in the addition of magnetic dipole–dipole interactions to the long-term balance, which delays the achievement of 〈M̃_eq〉. This is confirmed by the transient behavior of 〈M̃(t̃)〉 for systems with the clustering aggregation process (λ ≥ 15) observed in Fig. 13(a) and (b). However, this behavior shows the dependency of α and λ. In the same way, ϕ also affects them. This is because these parameters drive the clustering aggregation process. Therefore, the effect of each parameter on 〈M̃(t̃)〉 will be explained in terms of the corresponding clustering process, characterized by the microstructural properties previously analyzed.

Fig. 13 Comparison of temporal evolution of the magnetization, 〈M̃(t̃)〉, of a magnetic colloidal suspension with ϕ = 10⁻³, for λ = 15 (a) and λ = 75 (b) under different α values. The microstructure observed at t̃ = 1000 in both cases for α = 0.01, 1 and 1000 is shown below the magnetization curves.

Regardless of ϕ and λ, the transient behavior of 〈M̃(t̃)〉 is insignificant for very low and high α values. In the former case, such as α = 0.01, 〈N_c(t̃)〉 curves show chains growing as time goes on, but thermal fluctuations take them away from the rigid chain behavior, as shown by the 〈R_eff〉 curves, by forming rings (Fig. 8). The contribution of rings to 〈M̃(t̃)〉 is minimal because the sum of all its dipole moments approaches zero (negligible contribution). Additionally, with respect to flexible chains, the bending effect on them reduces the sum of their magnetic dipole moments in the magnetic field direction, causing a weak contribution. Thus, on average, 〈M̃(t̃)〉 ≈ 0.

For very high values of α, such as α = 1000, the transient behavior develops quickly and 〈M̃(t̃)〉 reaches its saturation value. The significant predominance of the magnetic field effect over thermal fluctuations yields a perfect alignment with it (see Fig. 8) and, consequently, a very rapid achievement of saturation magnetization. Thus, the clustering aggregation process does not greatly affect the transient behavior of 〈M̃(t̃)〉.

At intermediate α values, the transient behavior of 〈M̃(t̃)〉 is well-distinguishable, independent of ϕ and α. The ascending order of the 〈M̃(t̃)〉 curves with respect to α is explained by the stronger trend with higher α to align particles and clusters in the magnetic field direction. What needs to be explained is the mechanism that supports the transient behavior of 〈M̃(t̃)〉 due to the clustering aggregation process.

〈N_c(t̃)〉 and 〈n_c(t̃)〉 curves highlight the formation of short clusters in the nucleation stage. These behave like rigid chains, as shown by the 〈R_eff〉 curves (see Fig. 9). Thus, each short-dipolar chain can be modeled as a unique and stronger magnetic dipole moment. This encourages a greater magnetic field-induced torque on the chain, increasing its magnetic field alignment trend and augmenting 〈M̃(t̃)〉. However, the transient behavior of 〈n_c(t̃)〉 over time shows interaction between these chains, which deflects them from the aforementioned trend, resulting in the transition behavior of 〈M̃(t̃)〉.

Nevertheless, 〈N_c(t̃)〉 curves display the constant increase of the chains. The resulting structures, which continue to grow, gradually lose their rigid chain-like behavior, as observed in Fig. 9, and are distorted by thermal fluctuations (Fig. 8) for α = 1. These distortions develop differently along a chain, spending time to reach steady-state. Thus, this also delays the achievement of 〈M̃_eq〉.

This mechanism explains the origin of the transient behavior of 〈M̃(t̃)〉 for different α values due to the clustering aggregation process. However, the microstructure morphology can also influence 〈M̃(t̃)〉.

5.2.2 Microstructure morphology effect on the magnetization. With the progression from λ = 15 (Fig. 13(a)) to λ = 75 (Fig. 13(b)), 〈M̃(t̃)〉 increases faster. In order to understand this improvement, the observed microstructural differences must be analyzed. According to the 〈N_c(t̃)〉 curves (see Fig. 6), cluster sizes increase as λ increases. However, the most significant difference occurs for the low α regime, such as α = 1, as shown in Fig. 13.

For λ = 15, the enhancement of 〈M̃(t̃)〉 is mainly supported by flexible chains. According to the 〈R_eff〉 curves, their average behavior is close to the rigid chain behavior. This means that the contribution to 〈M̃(t̃)〉 is not negligible, especially given that they have great sizes, as shown by the 〈N_c(t̃)〉 curves. Furthermore, the presence of rings prevents further growth of magnetization.¹⁰⁰

For λ = 75, the chains have longer lengths. This complicates the formation of rings because the encounter of the ends of a chain via Brownian diffusion is more difficult. However, rings are present in the system, but in smaller quantities. Then, the small-ring population may lead to a higher contribution to 〈M̃(t̃)〉. In addition, chain flexibility allows for the formation of double chains, which consists of two steps: chain bending and zipper coupling. The last step ensures that the coupled chains have the same direction, which means that all their dipole moments strongly contribute to 〈M̃(t̃)〉 compared to the cases where double chains are not developed, such as λ = 15.

However, the formation of double chains is not necessary to enhance 〈M̃(t̃)〉. This is observed by analyzing the effect of ϕ on the transient behavior of 〈M̃(t̃)〉. For ϕ = 5 × 10⁻⁴ and high λ, the chains also become long and interact with each other, but do not form double chains. However, this also contributes to 〈M̃(t̃)〉. Their interaction forces them to have a similar orientation, improving 〈M̃(t̃)〉. In addition, the formation of short chains, as occurs for λ = 15, makes this effect difficult. Therefore, the enhancement of 〈M̃(t̃)〉 is not comparable to the case of high λ, similarly to the case of ϕ = 10⁻³.

Systems with ϕ = 10⁻⁴ exhibit the formation of short chains and short-lived rings in the low α regime. However, the contribution of this microstructure to 〈M̃(t̃)〉 can be substantial, depending on the value of λ.

5.2.3 Equilibrium magnetization. From the calculation of the equilibrium magnetization values 〈M̃_eq〉 using (30), saturation magnetization curves are constructed, showing the relationship between 〈M̃_eq〉 and the parameters α, λ, and ϕ. In Fig. 14, these curves are depicted and compared with the magnetization curve predicted by the Langevin function 〈M̃_eq〉 = L(α).

Fig. 14 Equilibrium magnetization curves of magnetic colloidal suspensions with ϕ = 10⁻⁴, ϕ = 5 × 10⁻⁴ and ϕ = 10⁻³. Each curve relates to a value of λ (λ = 15, 45, 75). The equilibrium magnetization 〈M_eq〉 is non-dimensionalized with the saturation magnetization ϕM_s of the suspension, 〈M̃_eq〉 = 〈M_eq〉/ϕM_s. The Langevin function, L(α), is used as a control curve.

In general, the influence of α on magnetization curves is well-described by the Langevin function. That is, the increase of α leads to an increase in 〈M̃_eq〉, reflecting the enhanced magnetic torque that promotes particle alignment along the magnetic field direction. This behavior is observed regardless of the values of λ and ϕ. In the absence of the clustering aggregation process, the Langevin function accurately predicts 〈M̃_eq〉, as occurs for λ = 5 and any value of ϕ (see Fig. 14). In contrast, the curves for λ > 5 show that the Langevin function underestimates 〈M̃_eq〉, mainly in the range 0.01 ≤ α ≤ 10, as observed in Fig. 14.

The influence of λ on the magnetization curves is similar to the former case. The increase in λ leads to an increase in 〈M̃_eq〉 (see Fig. 14). This is due to the formation of longer structures (see Fig. 6), as well as the lateral interchain interactions discussed previously.

In the same way, increasing the value of ϕ also increases the value of 〈M̃_eq〉. This occurs because a higher value of ϕ implies a reduction in the distance between the particles, which in turn promotes an enhanced clustering aggregation process. This is clearly evidenced in the magnetization curves for λ = 15, which exhibit higher values of 〈M̃_eq〉 as a result of the increasing value of ϕ (see Fig. 14).

In general, the influence of λ and ϕ manifests itself as deviations of the values of 〈M̃_eq〉 from those predicted by the Langevin function, as observed in Fig. 14. This suggests that the effects of both parameters may be combined into a single one. Several theoretical models developed to predict the equilibrium magnetization of a suspension adopt this approach by introducing the Langevin magnetic susceptibility, χ_L, defined as:

χ_L = 4λϕ. (31)

Therefore, magnetic colloidal suspensions that share the same value of χ_L should attain identical equilibrium magnetization, 〈M̃_eq〉, when α is held constant.

A subset of these models has been employed to compare their equilibrium magnetization predictions, denoted by M̃_eqm, with the 〈M̃_eq〉 values obtained in this work. The Langevin theory (LT) is the first comparative model. Although it does not depend on χ_L, this is a classical theory and serves as a valuable point of comparison. The second-order modified mean-field theory (MMF2) has been employed as the second comparative model.⁵⁷ This predicts equilibrium magnetization by defining an effective Langevin parameter α_e to be applied in the Langevin function, as follows:

M̃_eqm = L(α_e), (32)

As the third comparative model, the perturbation theory (PT) has been selected.⁵² By applying the non-dimensionalization process adopted in this work, the model yields the following equation for the equilibrium magnetization:

Since both models were developed to account for magnetic interactions between particles, they have successfully predicted equilibrium magnetization in concentrated suspensions characterized by low values of λ.

From Fig. 15, two important insights can be drawn. The first insight indicates that perturbation theory (PT) yields better predictive performance; however, it still deviates substantially from the 〈M̃_eq〉 values, especially within the range 0.01 < α < 10. As discussed previously, within this range, the formation of long chains and their lateral interactions occurs. This underscores the importance of thoroughly understanding these features to improve the accuracy of predictive models.

Fig. 15 Comparison between the values of the equilibrium magnetization, 〈M̃_eq〉, obtained via simulations with the predicted values, M̃_eqm, from theoretical models, such as the Langevin theory (LT), the perturbation theory (PT) and the second-order modified mean-field theory (MMF2) for the suspensions of λ = 75, ϕ = 10⁻⁴ and λ = 15, ϕ = 5 × 10⁻⁴.

The second insight relates to the values χ_L for the systems considered in Fig. 15. For both systems, χ_L = 0.03. Despite them sharing the same χ_L value, their corresponding PT curves in Fig. 15 differ significantly. This implies that the product of λ and ϕ does not result in a unified influence in diluted systems with moderate to strong magnetic interactions. This can be attributed to the large interparticle distances, which cause the relative alignment of the magnetic dipole moments to be governed predominantly by the strength of magnetic interactions, that is, by the λ value.

However, this does not mean that χ_L does not play a significant role in the characterization of dilute magnetic suspensions. To justify this, Table 2 shows the corresponding values of χ_L for any system simulated here. Then, the 〈M̃_eq〉 curves are plotted in Fig. 16, taking χ_L as constant. The corresponding curve for χ_L = 0.020 also represents the curves for the cases χ_L = 0.002 and χ_L = 0.010 (the Langevin function).

Table 2 Values of the Langevin magnetic susceptibility, χ_L, according to (31) and the values of the volume fraction, ϕ, and the dipolar coupling parameter, λ

χ L	ϕ	λ
0.002	10^−4	5
0.006	10^−4	15
0.010	5 × 10^−4	5
0.018	10^−4	45
0.020	10^−3	5
0.030 (1)	10^−4	75
0.030 (2)	5 × 10^−4	15
0.060	10^−3	15
0.090	5 × 10^−4	45
0.150	5 × 10^−4	75
0.180	10^−3	45
0.300	10^−3	75

---

Fig. 16 Equilibrium magnetization, 〈M̃_eq〉, curves as a function of the Langevin parameter, α, and characterized by their corresponding Langevin magnetic susceptibility, χ_L, values (according to Table 2).

Fig. 16 shows that the curves corresponding to χ_L < 0.090 do not exhibit a discernible pattern, whereas those with χ_L ≥ 0.090 tend to collapse onto a single curve, indicating a consistent behavior in this range. According to Table 2, this range corresponds to dilute suspensions with strong magnetic interparticle interactions. Hence, the range χ_L ≥ 0.090 will characterize the kind of suspension discussed in this work.

6 Conclusions

Using BD simulations, the influence of the clustering aggregation process on the time evolution of magnetization in magnetic colloidal suspensions is analyzed. These suspensions are characterized by different particle volume fractions (ϕ), applied uniform magnetic fields (α) and magnetic dipole–dipole interaction strengths (λ). A comprehensive description of the clustering aggregation process is attained using microstructural properties such as the nucleation-growth factor, 〈n_c(t̃)〉, mean cluster size, 〈N_c(t̃)〉, kinetic exponent, z, effective radius, 〈R_eff〉, and radial distribution function, g(r̃). Then, using the resulting description, analyses of the transient magnetization behavior, 〈M̃(t̃)〉, and equilibrium magnetization, 〈M̃_eq〉, are performed.

Simulations show a clustering aggregation process leading to a microstructure mainly composed of chains, which are flexible for α ≤ 1 (low α regime) and rigid for α > 1. However, the formation of other structures also takes place, for example, the formation of rings (low α < 1 regime) due to thermal fluctuations in the flexible chains, or the formation of double chains (α ≈ 1) due to the lateral bending of two interacting long flexible chains and their consequent zipper-like coupling.

The clustering aggregation process can be divided into a nucleation stage and a growth stage, both related to the initial particle–particle aggregation and the subsequent cluster–cluster aggregation, respectively. The behaviors of 〈n_c(t̃)〉 and 〈N_c(t̃)〉 over time reveal a higher aggregation rate as ϕ and λ increase. In contrast, the increase in α appears to have no significant effects on the clustering aggregation process. However, the case of α ≥ 1 leads to a delay in the aggregation process because ϕ and λ gradually take higher values. This delay is depicted by the decay of 〈N_c(t̃)〉. This suggests the presence of a certain effect that interferes with cluster aggregation. The analysis of microstructure morphology, 〈R_eff〉 and g(r̃) demonstrate that this is the result of lateral interaction between long rigid chains, which are paired for high ϕ and λ values.

Simulations with λ = 5 reveal a lack of clustering process regardless of the values of ϕ and α. Then, the Langevin function accurately predicts the 〈M̃_eq〉 curves for these cases. However, for higher λ values, the corresponding curves differ significantly from the previous case because the clustering aggregation process increases the magnetization of the suspension.

However, the clustering process also delays achieving 〈M̃_eq〉 when α has moderate values. According to 〈R_eff〉, the clustering process begins by producing short dipolar chains that on average behave like rigid chains. These can be modeled as short rigid chains with greater dipole moments that interact with each other. Then, as 〈N_c(t̃)〉 shows, these chains grow over time, further intensifying their dipole moments, and hence their interactions. These time-varying interactions delay the magnetization equilibrium state by producing a transient behavior. This behavior develops faster and has higher values as ϕ and λ increase.

As these short-dipolar chains grow, they gradually lose their rigid chain behavior. Then, the acquired flexibility allows the formation of rings (low α regime), which does not contribute to magnetization, or the formation of double chains, which hardly benefits magnetization because the dipole moments that compose them acquire the same orientation, which is nearby to the magnetic field direction (α ≈ 1). The presence of double chains is not necessary to enhance magnetization. If two long chains interact with each other, they will try to form a double chain. Despite the failure in this process, similar orientations acquired by them will be enough to increase the magnetization of the suspension. This is corroborated by comparing the equilibrium magnetization obtained from simulations with the theoretical predictions provided by the Langevin theory, perturbation theory, and second-order modified mean-field theory. The comparison reveals a pronounced discrepancy when the lateral interchain interactions develop, that is, in the range of 0.01 ≤ α ≤ 10. Finally, the Langevin magnetic susceptibility, χ_L, used in the comparative models serves as a parameter to characterize dilute suspensions with strong magnetic interparticle interactions. For χ_L ≥ 0.09, the nondimensionalized 〈M̃_eq〉 curves collapse onto a single one.

Author contributions

Luis R. Pérez-Marcos: conceptualization, formal analysis, methodology, writing – original draft, writing – review & editing. Ronal A. DeLaCruz-Araujo: funding acquisition, supervision, conceptualization. Heberth Diestra-Cruz: conceptualization, supervision, funding acquisition. Obidio Rubio: supervision, funding acquisition, resources. Ubaldo M. Córdova-Figueroa: resources, supervision. Glenn C. Vidal-Urquiza: supervision, formal analysis, methodology, funding acquisition, conceptualization, writing – review & editing.

Conflicts of interest

There are no conflicts of interest.

Data availability

The data for this article, specifically used for each figure, are available in the specific directories in Zenodo at https://doi.org/10.5281/zenodo.15671370.

Acknowledgements

This work was subsidized by CONCYTEC through the PROCIENCIA program within the framework of the contest “Basic Research Projects 2021-01”, according to contract no. 041-2021-FONDECYT. In addition, the authors acknowledge Boqueron HPCf from the University of Puerto Rico for the computational support provided for some simulations performed in this research work. Finally, Luis R. Pérez-Marcos thanks the NSF-EPSCoR Center for the Advancement of Wearable Technologies (CAWT) Graduate Fellowship for the partial financial support. NSF Grant No. OIA-1849243.

References

1. S. Odenbach, Colloidal magnetic fluids: basics, development and application of ferrofluids, Springer Science & Business Media, 2009, vol. 763.
2. R. E. Rosensweig, Ferrohydrodynamics, Courier Corporation, 2013.
3. J. Thomas, J. Appl. Phys., 1966, 37, 2914–2915.
4. R. Brown, A Brief Account of Microscopical Observations, etc., London(not published), 1827.
5. A. Einstein, Investigations on the Theory of the Brownian Movement, Courier Corporation, 1956.
6. W. Russel, Annu. Rev. Fluid Mech., 1981, 13, 425–455.
7. M. Matsuzaki, H. Kikura, J. Matsushita, M. Aritomi and H. Akatsuka, Sci. Technol. Adv. Mater., 2004, 5, 667–671.
8. M. Waarden, J. Colloid Sci., 1950, 5, 317–325.
9. E. Mackor, J. Colloid Sci., 1951, 6, 492–495.
10. R. Rosensweig, J. Nestor and R. Timmins, Materials associated with direct energy conversion, Institution of Chemical Engineers, 1965, pp. 104–118.
11. A. P. Philipse and D. Maas, Langmuir, 2002, 18, 9977–9984.
12. S. L. Tripp, S. V. Pusztay, A. E. Ribbe and A. Wei, J. Am. Chem. Soc., 2002, 124, 7914–7915.
13. K. Butter, P. Bomans, P. Frederik, G. Vroege and A. Philipse, Nat. Mater., 2003, 2, 88–91.
14. K. Butter, P. Bomans, P. Frederik, G. Vroege and A. Philipse, J. Phys.: Condens. Matter, 2003, 15, S1451.
15. M. Klokkenburg, C. Vonk, E. M. Claesson, J. D. Meeldijk, B. H. Erné and A. P. Philipse, J. Am. Chem. Soc., 2004, 126, 16706–16707.
16. J. Rabinow, Electr. Eng., 1948, 67, 1167.
17. K. D. Weiss, T. G. Duclos, J. D. Carlson, M. J. Chrzan and A. J. Margida, SAE Trans., 1993, 425–430.
18. J. D. Carlson, D. Catanzarite and K. St. Clair, Int. J. Mod. Phys. B, 1996, 10, 2857–2865.
19. B. F. Spencer Jr, G. Yang, J. D. Carlson and M. K. Sain, Proceedings of the second world conference on structural control, 1998, pp. 417–426.
20. P. Ilg, M. Kröger and S. Hess, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2005, 71, 051201.
21. N. K. Rai, G. Reddy, S. Ramanujam, V. Venkatraj and P. Agrawal, Def. Sci. J., 2009, 59, 239–251.
22. T. Sakurai and S. Morishita, Front. Mech. Eng., 2017, 12, 224–233.
23. F. Xiong, S. Huang and N. Gu, Drug Dev. Ind. Pharm., 2018, 44, 697–706.
24. P. Zhang, M. Kamezaki, K. Otsuki, Z. He, H. Sakamoto and S. Sugano, 2019 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM), 2019, pp. 400–405.
25. E. E. Etli and A. Akar, Medicine, 2022, 11, 934–941.
26. D. N. Rao Miriyala and P. Goyal, arXiv, 2022, preprint, arXiv:2201.09027 DOI:10.48550/arXiv.2201.09027.
27. L. R. Pérez-Marcos, R. A. DeLaCruz-Araujo, H. Diestra-Cruz, O. Rubio, U. M. Córdova-Figueroa and G. C. Vidal-Urquiza, Soft Matter, 2025 10.1039/D5SM00180C.
28. V. S. Mendelev and A. O. Ivanov, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2004, 70, 051502.
29. L. N. Donselaar, P. M. Frederik, P. Bomans, P. A. Buining, B. M. Humbel and A. P. Philipse, J. Magn. Magn. Mater., 1999, 201, 58–61.
30. S. Menear, A. Bradbury and R. Chantrell, J. Magn. Magn. Mater., 1984, 43, 166–176.
31. R. Chantrell, A. Bradbury, J. Popplewell and S. Charles, J. Phys. D: Appl. Phys., 1980, 13, L119.
32. G. Helgesen, A. Skjeltorp, P. M. Mors, R. Botet and R. Jullien, Phys. Rev. Lett., 1988, 61, 1736.
33. Z. Wang, C. Holm and H. W. Müller, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2002, 66, 021405.
34. A. Satoh, R. W. Chantrell, S.-I. Kamiyama and G. N. Coverdale, J. Colloid Interface Sci., 1996, 181, 422–428.
35. A. Satoh, R. W. Chantrell, S.-I. Kamiyama and G. N. Coverdale, J. Colloid Interface Sci., 1996, 178, 620–627.
36. M. Barrett, A. Deschner, J. P. Embs and M. C. Rheinstädter, Soft Matter, 2011, 7, 6678–6683.
37. P. De Gennes and P. Pincus, Phys. Kondens. Mater., 1970, 11, 189–198.
38. T. Tlusty and S. Safran, Science, 2000, 290, 1328–1331.
39. R. Chantrell, A. Bradbury, J. Popplewell and S. Charles, J. Appl. Phys., 1982, 53, 2742–2744.
40. C. F. Hayes, J. Colloid Interface Sci., 1975, 52, 239–243.
41. M. Klokkenburg, B. H. Erné, J. D. Meeldijk, A. Wiedenmann, F. A. V. Petukhov, R. P. Dullens and A. P. Philipse, Phys. Rev. Lett., 2006, 97, 185702.
42. M. Klokkenburg, B. H. Erné, V. Mendelev and A. Ivanov, J. Phys.: Condens. Matter, 2008, 20, 204113.
43. J.-J. Weis, Mol. Phys., 2005, 103, 7–10.
44. M. Ivey, J. Liu, Y. Zhu and S. Cutillas, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 63, 011403.
45. N. Mousavi, S. D. Khapli and S. Kumar, J. Appl. Phys., 2015, 117, 103907.
46. A. Darras, E. Opsomer, N. Vandewalle and G. Lumay, Eur. Phys. J. E: Soft Matter Biol. Phys., 2019, 42, 1–9.
47. A. Pshenichnikov and V. Mekhonoshin, J. Magn. Magn. Mater., 2000, 213, 357–369.
48. J. Li, Y. Huang, X. Liu, Y. Lin, L. Bai and Q. Li, Sci. Technol. Adv. Mater., 2007, 8, 448–454.
49. R. Kaiser and G. Miskolczy, J. Appl. Phys., 1970, 41, 1064–1072.
50. M. Shliomis, A. Pshenichnikov, K. Morozov and I. Y. Shurubor, J. Magn. Magn. Mater., 1990, 85, 40–46.
51. K. Morozov and A. Lebedev, J. Magn. Magn. Mater., 1990, 85, 51–53.
52. Y. A. Buyevich and A. Ivanov, Phys. A, 1992, 190, 276–294.
53. A. Pshenichnikov, J. Magn. Magn. Mater., 1995, 145, 319–326.
54. A. Y. Zubarev and L. Y. Iskakova, J. Exp. Theor. Phys., 1995, 80, 857–866.
55. A. Pshenichnikov, V. Mekhonoshin and A. Lebedev, J. Magn. Magn. Mater., 1996, 161, 94–102.
56. I. Abu-Aljarayesh and S. Migdadi, J. Magn. Magn. Mater., 1999, 191, 174–180.
57. A. O. Ivanov and O. B. Kuznetsova, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2001, 64, 041405.
58. S. Taketomi and R. D. Shull, J. Appl. Phys., 2002, 91, 8546–8548.
59. P. Langevin, J. Phys. Theor. Appl., 1905, 4, 678–693.
60. M. I. Shliomis, Soviet Physics Uspekhi, 1974, 17, 153.
61. M. Martsenyuk, Y. L. Raikher and M. Shliomis, Zh. Eksp. Teor. Fiz., 1973, 65, 834–841.
62. A. Sreekumari and P. Ilg, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2013, 88, 042315.
63. M. Fermigier and A. P. Gast, J. Colloid Interface Sci., 1992, 154, 522–539.
64. M. Hagenbüchle and J. Liu, Appl. Opt., 1997, 36, 7664–7671.
65. P. Domnguez-Garca, S. Melle, J. Pastor and M. Rubio, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2007, 76, 051403.
66. F. Martnez-Pedrero, M. Tirado-Miranda, A. Schmitt and J. Callejas-Fernández, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2007, 76, 011405.
67. A. O. Ivanov, Z. Wang and C. Holm, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2004, 69, 031206.
68. A. Satoh, Introduction to molecular-microsimulation for colloidal dispersions, Elsevier, 2003.
69. A. Satoh, Introduction to practice of molecular simulation: molecular dynamics, Monte Carlo, Brownian dynamics, Lattice Boltzmann and dissipative particle dynamics, Elsevier, 2010.
70. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, 2017.
71. W. Coffey and Y. P. Kalmykov, The Langevin equation: with applications to stochastic problems in physics, chemistry and electrical engineering, World Scientific, 2012, vol. 27.
72. J. C. Chen and A. S. Kim, Adv. Colloid Interface Sci., 2004, 112, 159–173.
73. A. Satoh, Modeling of magnetic particle suspensions for simulations, CRC Press, 2017.
74. G. G. Stokes, et al. , Math. Phys. Papers, 1851, 3, 1880–1905.
75. E. Guazzelli and J. F. Morris, A physical introduction to suspension dynamics, Cambridge University Press, 2011, vol. 45.
76. S. Chikazumi and C. D. Graham, Physics of ferromagnetism, Oxford University Press, 1997.
77. H. Nyquist, Phys. Rev., 1928, 32, 110.
78. R. Kubo, Science, 1986, 233, 330–334.
79. R. Kubo, Rep. Prog. Phys., 1966, 29, 255.
80. D. L. Ermak and J. A. McCammon, J. Chem. Phys., 1978, 69, 1352–1360.
81. E. Dickinson, S. A. Allison and J. A. McCammon, J. Chem. Soc., Faraday Trans. 2, 1985, 81, 591–601.
82. H. Goldstein, C. Poole and J. Safko, Classical mechanics, 2002.
83. D. J. Evans and S. Murad, Mol. Phys., 1977, 34, 327–331.
84. D. J. Evans, Mol. Phys., 1977, 34, 317–325.
85. G. I. Vega-Bellido, R. A. DeLaCruz-Araujo, I. Kretzschmar and U. M. Córdova-Figueroa, Soft Matter, 2019, 15, 4078–4086.
86. T. Vicsek and F. Family, Phys. Rev. Lett., 1984, 52, 1669.
87. S. Miyazima, P. Meakin and F. Family, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 1421.
88. J. H. Promislow, A. P. Gast and M. Fermigier, J. Chem. Phys., 1995, 102, 5492–5498.
89. G. Mériguet, M. Jardat and P. Turq, J. Chem. Phys., 2005, 123, 144915.
90. B. Akpinar, L. A. Fielding, V. J. Cunningham, Y. Ning, O. O. Mykhaylyk, P. W. Fowler and S. P. Armes, Macromolecules, 2016, 49, 5160–5171.
91. N. Kovalchuk and V. Starov, Adv. Colloid Interface Sci., 2012, 179, 99–106.
92. A. Tomilov, A. Videcoq, M. Cerbelaud, M. A. Piechowiak, T. Chartier, T. Ala-Nissila, D. Bochicchio and R. Ferrando, J. Phys. Chem. B, 2013, 117, 14509–14517.
93. F. L. Durhuus, L. H. Wandall, M. H. Boisen, M. Kure, M. Beleggia and C. Frandsen, Nanoscale, 2021, 13, 1970–1981.
94. J. Cernák, J. Magn. Magn. Mater., 1994, 132, 258–269.
95. A. Ben Salah, T. Ukai, L. Mingyuan, H. Morimoto and T. Maekawa, AIP Adv., 2023, 13, 055010.
96. D. Sohn, J. Magn. Magn. Mater., 1997, 173, 305–313.
97. J. Faraudo, J. S. Andreu, C. Calero and J. Camacho, Adv. Funct. Mater., 2016, 26, 3837–3858.
98. M. Mohebi, N. Jamasbi and J. Liu, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 5407.
99. J. S. Andreu, J. Camacho and J. Faraudo, Soft Matter, 2011, 7, 2336–2339.
100. L. Gutiérrez, L. De la Cueva, M. Moros, E. Mazaro, S. De Bernardo, J. M. De la Fuente, M. P. Morales and G. Salas, Nanotechnology, 2019, 30, 112001.

---

Footnote

† Electronic supplementary information (ESI) available: Videos of clustering aggregation for the suspensions with ϕ = 10⁻³ and λ = 75, for the case of α = 0.01 (Video 1), α = 1 (Video 2) and α = 1000 (Video 3). See DOI: https://doi.org/10.1039/d5sm00209e

---