On-the-fly coarse-graining methodology for the simulation of chain formation of superparamagnetic colloids in strong magnetic fields

The aim of this work is the description of the chain formation phenomena observed in colloidal suspensions of superparamagnetic nanoparticles under high magnetic fields. We propose a new methodology based on an on-the-fly Coarse-Grain (CG) model. Within this approach, the coarse grain objects of the simulation are not fixed a priori at the beginning of the simulation but rather redefined on the fly. The motion of the CG objects (single particles or aggregates) is described by an anisotropic diffusion model and the magnetic dipole-dipole interaction is replaced by an effective short range interaction between CG objects. The new methodology correctly reproduces previous results from detailed Langevin Dynamics simulations of dispersions of superparamagnetic colloids under strong fields whilst requiring an amount of CPU time orders of magnitude smaller. This substantial improvement in the computational requirements allows the simulation of problems in which the relevant phenomena extends to time scales inaccessible with previous simulation techniques. A relevant example is the waiting time dependence of the relaxation time T_2 of water protons observed in Magnetic Resonance experiments containing dispersions of superparamagnetic colloids, which is correctly predicted by our simulations. Future applications may include other popular real-world applications of superparamagnetic colloids such as the magnetophoretic separation processes.

The aim of this work is the description of the chain formation phenomena observed in colloidal suspensions of superparamagnetic nanoparticles under high magnetic fields. We propose a new methodology based on an on-the-fly Coarse-Grain (CG) model. Within this approach, the coarse grain objects of the simulation and their dynamic behavior are not fixed a priori at the beginning of the simulation but rather redefined on-the-fly. The motion of the CG objects (single particles or aggregates) is described by an anisotropic diffusion model and the magnetic dipole-dipole interaction is replaced by an effective short range interaction between CG objects. The new methodology correctly reproduces previous results from detailed Langevin Dynamics simulations of dispersions of superparamagnetic colloids under strong fields whilst requiring an amount of CPU time orders of magnitude smaller. This substantial improvement in the computational requirements allows the simulation of problems in which the relevant phenomena extends to time scales inaccessible with previous simulation techniques. A relevant example is the waiting time dependence of the relaxation time T2 of water protons observed in magnetic resonance experiments containing dispersions of superparamagnetic colloids, which is correctly predicted by our simulations. Future applications may include other popular real-world applications of superparamagnetic colloids such as the magnetophoretic separation processes.

I. INTRODUCTION
In recent years, work in coarse grain models for the description of soft matter and biomolecular systems is experiencing a remarkable outburst [1]. The reason is that the description of these systems at experimentally relevant time and length scales requires inclusion of phenomena occurring at very different scales. The objective of coarse grain (CG) models is thus to retain sufficient molecular or nanoscale detail and yet remain amenable of simulation up to macroscopic time scales. Many approaches have been developed to construct CG models of different kinds of soft matter systems. For example, in the case of polymers, there is a long tradition of using CG models and the field is sufficiently mature so that there are systematic and rigorous approaches to build up CG models from accurate atomistic descriptions [2]. Also, in the field of biomolecular simulations, there are important developments such as the MARTINI force field [3] which allow the simulation of difficult problems such as the behavior of lipid vesicles or protein folding at millisecond or even larger time scales. New advances include also simulation packages specially designed for CG models of soft matter such as ESPResSo [4].
Our interest here is the development of an improved CG model for a specific problem which is still difficult to simulate, namely the assembly of superparamagnetic colloids under strong magnetic fields. Superparamag- * Electronic address: jandreu@icmab.es netic colloids are made of small nanoparticles of magnetic material (typically 5-10 nm iron oxide nanocrystals) embedded in a nonmagnetic matrix (typically polymers or silica) [5]. These particles have no magnetic dipole in absence of magnetic field but they develop very high magnetizations in the presence of a magnetic field, similar to those obtained with ferromagnetic materials. This highly tunable response and the possibility to functionalize their surface make these materials very interesting for applications such as capture and removal of biomolecules and pollutants, NMR contrast agents, and many others [6][7][8].
Our work is motivated by the difficulties encountered in modeling different kinds of real experimental situations involving superparamagnetic colloids. A relevant example is provided by the experiments by Chen et al. [9] of a dispersion containing superparamagnetic colloids designed as contrast agents for magnetic resonance imaging (MRI). In these experiments, a strong uniform magnetic field was applied to the dispersion. It was observed that the transverse relaxation time T 2 of protons in water changed with time, an effect which was attributed to the formation of chains of superparamagnetic colloids. In fact, the kinetics of chain formation was estimated from these experiments, spanning time scales from 10 to 10 3 seconds or more. Another relevant example is magnetophoresis [10], which is the motion of magnetic particles under an external magnetic gradient. Experimental evidence shows that the formation of chains induced by the external field speeds up the magnetophoresis process [11], which is orders of magnitude faster than that observed in absence of chain formation [12]. It is worth noting that in these experiments, chains dissolve almost immediately after removal of the external field, as should be expected since superparamagnetic nanoparticles have no dipole in absence of magnetic field. In this respect, these systems are very different from the widely studied (and simulated) dispersions of dipolar particles, which are able to form structures in the absence of an external magnetic field due to the interaction of their permanent dipoles [13][14][15].
The standard approach for simulation of chaining processes in magnetic colloids is the use of Langevin Dynamics (LD) simulations (see for example [16,18]). This technique allows the inclusion of particle-particle interactions, thermal noise and the friction due to the fluid. The resolution of the simulation technique is typically in the nanoseconds scale. Simulation runs up to a few seconds are possible, but they are highly intensive, requiring the use of parallel computing during several weeks [18]. These CPU requirements make this simulation technique unsuitable for the study of the problems mentioned above.
The need to account for microscopic time and length scales but also reach macroscopic time scales at low computational cost has motivated us the development of a new simulation strategy based on an on-the-fly CG procedure. The methodology which will be developed in this paper is a generalization of two procedures proposed in previous works: the method proposed by Miguel and Satorras [19] to study aggregation processes and the method proposed by Schaller et al. [20] to study magnetophoresis.
In the methodology proposed here, one starts by simulating the motion and interaction between individual colloids. As the simulation advances, colloids form chains due to the magnetic dipole-dipole attraction induced by the high external magnetic field. The motion of each particle inside a chain is not simulated explicitly. In our methodology, these chains are considered individual coarse grain (CG) objects which move following certain effective rules and interact (and possibly aggregate) with other CG objects or single individual particles. In this way, the CG objects of the simulation are not fixed a priori at the beginning of the simulation but rather redefined on-the-fly. Thus, we adjust the resolution of the calculations during the simulation run, allowing for the possibility of much longer simulation runs requiring less computer power. Preliminary simulation results and comparison with experiments, presented in a previous work [9], demonstrated the feasibility and utility of our novel approach. Here we will discuss in detail the physical basis of the model, the simulation methodology and detailed comparison with more standard Langevin simulation techniques. All simulations of our model were performed employing the MagChain program, a C++ application developed in-house, which is freely available for use of researchers. The code, its documentation and usage examples can be found available for download at our web page [21].
The paper is organized as follows. In Section II, we describe the modeling of the system under study and the simulation technique. In Section III we validate the new methodology by (i) comparing our results with those obtained employing standard LD simulations (ii) discussing the effect of choosing other approximations for the diffusion coefficients and the effective interaction of the CG objects and (iii) discussing the applicability of our methodology comparing with experimental results. The conclusions are presented in section IV and some technical issues are detailed in the Appendix.

II. FORMULATION OF THE MODEL
The system which we are interested to describe is a colloidal dispersion of N p superparamagnetic spherical particles of diameter d in a volume V and volume fraction φ 0 : In absence of external magnetic field, the particles have no magnetic dipole and there is no formation of chains (no aggregation induced by the magnetic field). In the presence of a magnetic field H, the superparamagnetic particles acquire a certain magnetization M (H). Since we are particularly interested in the case of very strong magnetic fields (as in the experiments of Ref. [9] for example), we consider that the particles have a magnetic dipole moment m s (corresponding to the saturation magnetization M s of the particles) pointing in the direction of the applied magnetic field (which we will take as the z axis). The strength of the magnetic interaction between particles as compared with thermal energy can be characterized by the magnetic coupling parameter Γ defined as: The behavior of superparamagnetic colloids under external fields is controlled by the values of these two parameters (φ 0 and Γ). In this paper, we are interested in situations (values of φ 0 and Γ) in which the external field induces formation of linear chains of colloids which grow irreversibly with time. Irreversible growth of linear chains has been found in simulations and experiments investigating the ranges of Γ between 40 and 3×10 3 and φ 0 < 0.15 [16][17][18]. However, other structures are found at different ranges of φ 0 and Γ. For lower values of Γ, an equilibrium state is possible, in which colloids aggregate in linear (non branched) chains with an equilibrium length given by φ 0 e Γ−1 [18]. In the opposite situation of larger values of φ 0 and Γ, different aggregate structures can be found, including thick chains (obtained from lateral aggregation of linear chains), bundles and more complex fibrous structures [17]. All these more complex situations, different from irreversible growth of linear chains, will be left for future extensions of the model.
As key ingredients to retain the underlying physics of irreversible chain growth, we consider both the diffusive motion of particles and chains and their respective magnetic and steric interactions. The main approximations of our model will be to ignore the details of the particles forming the chains and to replace the actual magnetic dipole-dipole interaction by an effective short-range interaction, less demanding from the computational point of view.
Our model to study the kinetics of chain formation in these systems consists of CG objects which are chains made of s particles, including the case s = 1 which corresponds to a single particle. The first ingredient of the model is the diffusion coefficient of the CG objects. For single, isolated particles (s = 1) we have: where η is the viscosity of the fluid. A chain containing s > 1 particles exhibits anisotropic diffusion, characterized by a diffusion coefficient D s in the direction parallel to the long axis of the chain and D ⊥ s in the directions perpendicular to the long axis. There are several possibilities for the analytical form of these diffusion coefficients, depending on the exact geometry assumed for the chains and the degree of approximation of the calculation. Here, in order to keep the model as simple as possible, we consider the following expressions valid for elongated objects (slender body theory [22]): Strictly speaking, Eqs. (4) and (5) are valid only for large s. Therefore, we employ Eqs. (4) and (5) for chains with s > 2 and use a simple interpolation between the diffusion coefficients corresponding to CG objects with s = 1 (Eq. 3) and s = 3 (Eqs. (4) and (5)) for chains with s = 2. Such a choice gives results indistinguishable from those provided by more sophisticated and accurate expressions of the diffusion coefficients (see Section IIIB). The second ingredient of the model is the definition of the effective interaction between CG objects. A CG object of length s interacts with other CG objects through an excluded volume interaction (hard core) corresponding to a cylinder of length s × d and diameter d. They also interact through dipole-dipole interactions. In order to simplify and speed up the simulations, we have replaced the actual dipole-dipole magnetic interaction between colloids by an effective, short range interaction between the CG objects. This interaction is defined as follows. For a given CG object, we define two spherical attractive regions of radius r a (s) (which depend on the length of the chain s) located at the two ends of the chain. As illustrated in Fig. 1, these regions are designed to mimic the region at which the magnetic attraction between a chain of particles (magnetized in the z direction) and an incoming test dipolar particle is equal or stronger than the thermal energy k B T . The values of r a (s) are calculated by finding the distance in the z axis at which the magnetic interaction energy E mag between a chain of s particles and a single test particle is equal to −k B T . Therefore, r a (s) is given by the solution of: The results of Eq. (6) for different values of Γ are also shown in Fig. 1. Once we have defined the range of the interaction, we need to define the strength of this interaction. In order to keep our model as simple as possible, we simply assume that all events in which a CG object enters into the interaction region of another CG object will lead to instantaneous aggregation. This rule has been employed previously in the interpretation of experimental results and it has been suggested by direct observations of chain formation under a microscope (see Refs. [23,24]). As we will show in the next section, this rule reproduces correctly the results obtained from LD simulations in which the magnetic interaction is computed accurately. The sensitivity of the results to the choice of r a will be also discussed in Section IIIB.
Once the basic ingredients for the model (rules for motion and interaction) are defined, it is necessary to specify the algorithm for the numerical solution of the model. In our case, the diffusive motion of the CG objects is simulated using the brownian dynamics technique [25]. At each time step ∆t a random displacement in each direction is generated with a gaussian distribution with zero mean and variance 2D s ∆t, where D s is the diffusion coefficient of the object (single particle or chain) in the direction of motion (x, y or z). Also, at each time step the distances between CG objects are checked in order to detect penetration of a CG object inside the region of aggregation of another CG object, as explained above, or to detect possible overlaps between them. In the case of aggregation of two CG objects, a new CG object is created (and the two previous CG objects are erased from the simulation) with length s, obtained from adding the lengths of the two aggregating chains and located at the center of mass of the aggregating CG objects. In the case of overlap between two CG objects without penetration into the aggregation region, we consider that the two chains collide. In this case, the moving CG object is placed in contact with the other one (without overlapping) at the collision coordinates defined by the trajectory previously followed (see Fig. 2). Finally, it should be noted here that the selection of an appropriate time step ∆t for the simulation is a crucial issue. A detailed discussion on the selection of ∆t is given in the Appendix.
Hence, a typical simulation run is as follows. The simulation starts from a pre-equilibrated system containing N p colloids (CG objects with s = 1). As the simulation goes on, colloids aggregate and chains with increasing values of s appear. Consequently, the number of CG objects of the simulation decreases with time and the simulation speeds up as the time advances, as we will discuss in de-tail in the following section. As a simulation output, we obtain the number of chains containing s colloidal particles at time t, n s (t). During the simulation, we also monitor the time evolution of the average number of colloidal particles in a chain N (t) defined as in [16,19]: and the probability of finding an aggregate of size s at a given time, defined as:

A. Comparison with Langevin Dynamics simulations
Our objective in this section is to compare the performance and results obtained using the model described in Section II with standard LD simulations of the same system. Briefly stated, Langevin Dynamics simulations consist on solving the Newton equations of motion for each particle taking into account external forces, the interaction forces between particles (magnetic and steric), the viscous drag from the solvent and a stochastic force arising from the thermal noise due to the fact that the  (1) and (2), ρp is the density of a single colloid, d is its diameter and D1 its diffusion coefficient. T is the temperature of the dispersion and η the viscosity of the solvent (water) at this temperature.   system is at a given temperature T . This comparison between our simplified CG method and more detailed LD simulations will help to clarify the validity of the approximations introduced in our model, as described in the previous section. In order to perform a significant comparison between the new procedure and the standard LD simulation technique, we have selected two cases with very different magnetic coupling Γ which were studied in previous works. The details for these systems are sum-marized in Table I and were denoted as Case 1 and Case 2.
Let us consider first Case 1, which corresponds to a dispersion of 100 nm superparamagnetic colloids at a volume fraction φ 0 = 5.23 × 10 −4 which have a magnetic coupling parameter Γ = 40 at saturation (i.e. under strong magnetic fields). This system was studied employing LD simulations in Ref. [18] by using the standard LAMMPS simulation package [27] (version 21May2008).  Now, we will compare these published results obtained with the standard LD technique with our CG methodology described in the previous section. These simulations will be denoted as LD40 (the Langevin Dynamics case) and CG40 (our Coarse Grain methodology). The parameters employed in the numerical algorithm are given in Table II. For completeness, we also give the parameters employed in the LD simulations. It should be noted that the LD simulations require a small time step (of the order of the ns). This small time step is needed in order to avoid a simulation crash in simulations involving chains of colloids, since the motion of colloids inside a chain involves very small displacements which need to be resolved with high precision. In our CG methodology (which do not consider the structure of chains), we can use much larger time steps as shown in Table II. A detailed discussion on the selection of the time step in our methodology is given in the Appendix. Typical snapshots of the simulation are shown in Fig. 3 and the results for N (t) are shown in Fig. 4. The snapshots illustrate the different resolution employed in the CG40 and LD40 simulations. As seen in the snapshots, the LD40 simulation resolves the individual particles making up the chains whereas the chains are structureless in the CG40 simulation. It should be noted that the chains obtained in the LD40 simulation are almost perfectly linear and are not significantly different from the coarse-grain objects of the CG40 simulation. As shown in Fig.4, the values of N (t) obtained from both simulations (CG40 and LD40) are in excellent agreement. For example, at t = 1s, the mean aggregate size for the CG40 simulation is N CG =12.10 and the value calculated from LD simulations is almost identical, N LD =12.14. Therefore, we can conclude that the simplifying approximations included in the new methodology (particularly those regarding the calculations of particleparticle magnetic interaction) do not affect the average size of chains.
We have performed a more detailed comparison between both approaches by comparing the distribution of chains of size s at certain times. In Fig. 5 we compare the corresponding probability distribution (defined in Eq. (8)) at t =1s obtained from LD40 and CG40 simulations. The agreement between both results after 1s is remarkable, and only slight differences are observed. As shown in Fig. 5, the distribution of chain sizes is very broad, with significant probabilities of finding chains well above and well below the average length (including isolated particles). Now, let us consider the case denoted as Case 2 in Table I corresponding to one of the samples considered in the experiments in [9]. In this case, the particles have larger saturation magnetization (Γ = 247) but the dispersion is more diluted (φ 0 = 4.64 × 10 −6 ). We have performed Langevin Dynamics simulations as well as simulations employing our new methodology, as in Case 1. These simulations will be denoted as LD247 and CG247, respectively. A list of relevant simulation parameters for LD247 and CG247 simulations is given in Table II. The results obtained for N (t) are also given in Fig. 4. The results for the probability distribution p(s; t) at t = 5s are shown in Fig. 6. Again, we obtain a good agreement between the predictions of both simulation methodologies, both in the average size of chains and in the probability distribution of chain sizes. Table II, both for the case with Γ = 40 and Γ = 247, the computational cost of the CG simulation technique is much lower than the corresponding LD simulation. For example, we note that a production run of 6 seconds for the LD247 simulation requires about 4.5 days of calculations, with the program running in parallel in a 8-Core AMD Opteron Magny Cours 6136 processor. In contrast, the CG247 simulation requires less than 4 hours to simulate the same physical time using a single core of the same processor employed in the LD247 run. In addition, we can reach surprisingly long time scales in our CG247 simulation with a very low computational cost (see Table II). In simulation CG247, we reach simulated times up to 10 3 s in a 1 day calculation, a time scale two orders of magnitude larger than that accessible using Langevin Dynamics simulations.

As shown in
The CPU costs shown in Table II demonstrate that the new simulation technique allows us to perform simulations of the two systems considered here with an extremely reduced computational effort as compared to Langevin Dynamics simulations. Moreover, it is also important to notice that the required CPU time for the CG approach to simulate a certain time interval is reduced during a simulation, since the number of CG objects decreases as the simulation advances. This effect is clearly shown in Fig. 7 for the CG40 simulation. We also show that the rate between the elapsed CPU time and the corresponding real time simulated depends linearly with the number of CG objects (see inset Fig. 7). It should be noted that in the LD simulations the opposite effect was observed; the fact that the individual motion of the particles inside the chains are fully resolved makes LD simulations increasingly inefficient as time goes on. B. Further discussion on the approximations of the model

Diffusion Model
As it has been already mentioned, one of the two key ingredients of the CG approach is the diffusion model adopted to describe the motion of the coarse grain objects. In order to check the possible influence of the model selected, we have also performed simulations with a different diffusion model proposed by Tirado et al. [28] in which they describe the translational motion of right circular cylinders also accounting for the so called endeffects. Following the same approach than in [23], we have used the expressions: where γ and γ ⊥ are the end-effect functions defined as: γ ⊥ (s) = 0.84 + 0.18 We have computed the mean aggregate size N using  4) and (5) and represented by solid circles) and the cylinder approximation with end-effects (corresponding to Eqs. (9) and (10) and represented by open circles) when computing the mean chain size of the aggregate N for the CG40 system.
both diffusion models for the CG40 system and the results obtained are plotted in Fig. 8. As it can be seen from these results, no significant differences are found in the average number of particles N obtained with both diffusion models. For this reason we can conclude that both models are suitable for the description of the diffusive motion of the chain-like aggregates in such systems and our selection of the elongated rod model (Eqs. (4) and (5)) instead of Eqs.(9)-(12) for the simulations was based on its major simplicity.

Effective interaction. Attraction radius
As explained in detail in Section II, we have defined the aggregation regions for each CG object as the surrounding space in which the magnetic interaction energy between the CG object and a dummy single particle is equal or smaller than −k B T . As shown in Fig. 1, the attraction radius r a defining this region depends on the size of the considered aggregate and on the magnetic coupling parameter Γ (see Eq. 6). It is observed that for small chains the attraction radius increases with their size and tends to a constant value for larger chains (the addition of a new particle into the same aggregate does not significantly contribute to the interaction magnetic energy). Here, we would like to demonstrate the importance of accounting for the s dependence of r a (s) in the simulations.
To this end, we have performed two additional simulations in which the s dependence of r a (s) is ignored. A first simulation (denoted as CG40-min) corresponds to a repetition of the simulation CG40 of the previous section (see Table II) but using r a = 1.46d for all chains. We also performed another simulation (denoted as CG40-max) in which we employed the value r a = 2.20d for all chains. These values correspond to the minimum and maximum values of r a (s) employed in the original CG40 simulation (see Fig. 1).
All these approaches give us different dynamics of the system as it is shown in Fig. 9 where the mean aggregate size N is plotted as a function of time together with the corresponding LD results. We observe that the CG40 simulation evolves from an initial behavior close to the CG-min simulation to a behavior closer to the CGmax simulation. Analogous calculations for the CG247 system (not shown here) exhibit identical behavior. In consequence, is important to take into account the full r a (s) dependence in the simulations as described in our formulation of the model in Section II. C. An example of practical application: Chain growth and T2 measurements Our objective in this subsection is to illustrate the applicability of the methodology developed here in situations of interest for applications of superparamagnetic III: Set of parameters used for numerical integration in the Coarse Grain simulations of the same colloids described in Case 2 in Table I but for different volume fractions φ0. Np is the number of particles in the simulation, Lz and Lx = Ly are the sizes of the simulation box (in units of particle diameters) in the directions parallel and perpendicular to the magnetic field, respectively. ∆t is the time step and t f is the total simulated time. As in Table II, Table III. colloids. As an example, let us consider the use of superparamagnetic colloids as contrasts agents in magnetic resonance imaging. An important issue in this application is the possibility of chain formation of colloids due to the strong magnetic fields applied in the experiments. The formation of chains of colloids in the sample increases the transversal relaxation time T 2 of protons, which is an undesired effect in practice. In a previous work, we have employed a preliminary version of our simulation code to analyze this possibility in MRI [9]. We have found that under conditions of interest for MRI, significant chaining occurs. We would like to discuss here the results of our simulations as well as compare our results with the experiments in a more direct way than the preliminary simulations presented in Ref. [9]. The system considered here is a dispersion of superparamagnetic colloids in water with the physical properties of Case 2 in Table I but now we have considered 4 different values of the initial volume fraction of colloids, according to the experiments in Ref. [9]. The corresponding volume fractions are given in Table III. The simulations of these systems, performed with the methodology discussed in section II have been labeled as CG247-00, CG247-01, CG247-02 and CG247-03, respectively and all the technical details are given in Table III (note that the CG247-02 simulation in this table is identical to the simulation CG247 of Table II).
The results for the average number of particles in a chain N (t) are given in Fig. 10 for the time scales relevant in the experiments of Ref. [9]. In all cases we observe significant chain formation even in the case of the smallest concentration. Of course, the kinetics of chain formation is observed to slow down as the concentration of colloids decreases.
The formation of chains has direct impact on the transversal relaxation rate 1/T 2 of water protons. Initially (t = 0), the relaxation rate of water protons 1/T (0) 2 is determined by the presence of a random distribution of isolated (dispersed) colloids. As time goes on, chains form and modify the T 2 response of the surrounding water protons. Therefore, the experimentally measured T 2 at a given instant t depends on the distribution of chain sizes at that time t. As proposed in Ref. [9], we can give a theoretical prediction of the relaxation rate 1/T 2 (t) from our simulation results by computing the following average: In Eq. (13), 1/T (s) 2 is the relaxation rate of water protons near a colloid forming part of a chain containing exactly s colloids and n s (t) is the number of chains of size s at time t, as defined in Section II. Our simulation results provide n s (t) whereas the calculation of 1/T were given in Fig. 9 of Ref. [9]. These results can be well fitted to an analytical expression of the form: where a fit to the calculations in Ref. [9] gives a = 0.0415 and b = 0.45. Now, making use of such a fit in Eq. (13) and the values of n s (t) obtained from simulations CG247-00, CG247-01, CG247-02 and CG247-03, we can make a theoretical prediction for the relaxation rate 1/T 2 (t).
The results are compared in Fig. 11 with experimental results. The simulations show a remarkable agreement between theory and experiments for times corresponding to mean chain length N larger than 50 colloids. It should be noted that in the case of very long chains, the measurements are not reliable due to sedimentation effects [9].  Table III and Eqs. (13) and (14). Symbols correspond to experimental data extracted from Fig.5a in Ref [9].

IV. CONCLUSIONS
In this paper, we have presented a new on-the-fly coarse grain model to describe the chaining phenomena observed in dispersions of superparamagnetic colloids under strong external magnetic fields. We report simulation results with the new methodology, which show good agreement with those obtained from more detailed Langevin Dynamics (LD) simulations. The great advantage of the new methodology presented here is its low computational cost in terms of CPU time. As a consequence, we are able to run longer simulations, reaching time scales not accessible in LD simulations. In order to illustrate the applicability of the code in experimentally relevant situations, we have considered the waiting time dependence of the relaxation rate 1/T 2 of water protons observed in magnetic resonance experiments of dispersions of superparamagnetic colloids [9]. Experimental results corresponding to waiting times from 1 to 10 3 s were correctly predicted by our simulations.
The model, in its present formulation, cannot be applied to situations more complex than irreversible chain growth. However, it seems possible to expand the model to consider other situations of interest. A first generalization could involve the inclusion of lateral interactions between the chains [29], which are responsible for the formation of thick chains, observed at volume fractions larger than those considered here [17]. Optical microscopy observations [11] also show the formation of thick chains and bundles in magnetophoresis experiments (motion of magnetic particles under magnetic gradients). Hence, the inclusion of lateral interactions and deterministic motion of the CG objects will be needed in order to extend our model to study magnetophoresis. Another interesting extension, which is now under way, is the inclusion of the possibility of chain breaking due to thermal fluctuations, a mechanism relevant at low values of the magnetic coupling parameter Γ. This extension of the model will allow us to study in depth the equilibrium state described in Refs. [18,30].
A final improvement to the model could be taking into account the full magnetic response M (H) of the particles in the simulation, in order to simulate situations in which the external magnetic fields are not strong enough to saturate the magnetic colloids. This is a typical situation in many published experimental studies of aggregation of magnetic colloids (see for example [23,24]), which focus on the linear magnetic response regime of the colloids.
where we have defined the characteristic diffusion time τ ≡ d 2 /D. In general, one selects a timestep ∆t which results in a displacement smaller than the relevant length scales of the problem (typical separations between particles, range of interaction forces,...). In our model, the length scale of interactions is given by the radius of the attraction zones (see Fig.1). The typical diffusive displacement corresponding to the selected ∆t (Eq.(A1)) has to be smaller than the radius of the attraction zone r a of the CG objects. In this way, CG objects will correctly explore the attraction zone of other surrounding CG objects. As it is shown in Fig.1c, the size r a of the attraction zone depends on the chain length s (the smallest value of r a corresponds to s = 1) and on the coupling parameter Γ. The dependence on Γ is strong, so one has to take into account this fact in selecting ∆t.
In order to check the effect of ∆t in the results of our simulations at Γ = 40 and Γ = 247, we have repeated the CG40 and CG247 simulations with three different timesteps: ∆t = 1.00τ , 0.10τ and 0.01τ . These timesteps ∆t correspond to typical displacements of 2.4d, 0.8d and 0.24d respectively (see Eq.(A1)). The results of these simulations for the average number of particles in a chain are shown in Fig.12. As it can be observed, the effects of the selection of the timestep are critical for the CG40 system and irrelevant for the CG247 system. In order to understand the effect of these different ∆t , one has to compare the obtained for each ∆t with the values of r a (s = 1) calculated for Γ = 40 and Γ = 247 (see Figure 1c). In the CG40 simulation, the smallest radius of a CG attraction zone is r a (s = 1) = 1.46d (see Figure  1c, case Γ = 40), so the attraction sphere has a diameter similar to the displacement = 2.4d obtained for ∆t = 1.0τ . In this case, colloids cannot explore properly the attraction zones and many chain formation events are lost during the simulation run. The other two ∆t give almost identical results since in both cases is smaller than r a (s = 1), and attraction zones are explored properly. In the CG247 simulations, we observe almost identical results for the three selected ∆t. In this case, we have r a (s = 1) = 2.89d (see Figure 1c, case Γ = 247), which is larger than the corresponding typical displacements . As a general rule, if the timestep selected is too large, the chains cannot properly explore the binding sites of the surrounding chains and the chaining process is not correctly simulated.