Multicore-based ferrofluids in zero field: initial magnetic susceptibility and self-assembly mechanisms

Abstract

The necessity to improve magnetic building blocks in magnetic nano-structured soft materials stems from a fascinating potential these materials have in bio-medical applications and nanofluidics. Along with practical reasons, the interplay of magnetic and steric interactions on one hand, and entropy, on the other, makes magnetic soft matter fundamentally challenging. Recently, in order to tailor magnetic response of the magnetic particle suspensions, the idea arose to replace standard single-core nanoparticles with nano-sized clusters of single-domain nanoparticles (grains) rigidly bound together by solid polymer matrix – multicore magnetic nanoparticles (MMNPs). To pursue this idea, a profound understanding of the MMNP interactions and self-assembly is required. In this work we present a computational study of the MMNP suspensions and elucidate their self-assembly and magnetic susceptibility. We show that depending on the magnetic moment of individual grains the suspensions exhibit qualitatively distinct regimes. Firstly, if the grains are moderately interacting, they contribute to a significant decrease of the remanent magnetisation of MMNPs and as such to a decrease of the magnetic susceptibility, this way confirming previous findings. If the grains are strongly interacting, instead, they serve as anchor points and support formation of grain clusters that span through several MMNPs, leading to MMNP cluster formation and a drastic increase of the initial magnetic response. Both the topology of the clusters and their size distribution in MMNP suspensions is found to be notably different from those formed in conventional magnetic fluids or magnetorheological suspensions.

---

1 Introduction

Nowadays, the idea of using smart nano-textured materials that we can thoroughly control and manipulate with external stimuli forms the basics for novel medical applications. In particular, when it comes to magnetic fields that are not interfering with any biological or physiological processes,¹ magnetic soft matter systems have gained a lot of attention for their potential in drug targeting^2–8 and magnetic hyperthermia.^9–12 The beginning of the magnetic soft matter field can be associated with the synthesis of ferrofluids back in the middle of the twentieth century.¹³ Today they are regarded as, probably, the simplest, albeit not fully understood and exploited, example of magneto-controllable material. Ferrofluids (or magnetic fluids) consist of two components: magnetic single-domain nanoparticles suspended in liquid magneto-passive carriers. It is the polydispersity of magnetic nanoparticles, peculiarities of their colloidal stabilisation and the resulting span of time-, length- and energy-scales that make ferrofluids difficult to manipulate in a highly controlled manner. In fact, any structures that nanoparticles form in ferrofluids are very sensitive to thermal noise and mechanical perturbations. Instead, it was recently suggested to replace single-domain magnetic particles by preassembled complex objects such as supramolecular polymer-like structures made of magnetic monomers,^14–24 magnetic nano- and micro-gels,^4,25,26 or magnetic multicore nanoparticles (MMNPs).^27–30 The suspensions of the latter were recently even addressed as “bio-ferrofluids”.³¹ This name is not just metaphoric, as MMNPs already gave rise to innovations in drug delivery,^32,33 magnetic particle imaging,^27,34 magnetic hyperthermia cancer treatment^35–38 and immunoassays.^39,40 Along with biomedical applications MMNPs have been used in nanorheology.^41–46 In all these works, one can think of a single MMNP as of a cluster of single-domain magnetic nanocrystals (below in this work addressed as “grains”) embedded in a polymer or other non-magnetic rigid matrix.^47,48 While the grains typically have a characteristic linear size of the order of 10 nm, the size of MMNPs can range from tens to a few hundred nanometers.

Despite such a wide range of promising applications, the underlying mechanisms, intrinsic and field-induced interactions of the MMNPs in suspensions are far from being understood. There are several theoretical works elucidating the impact of magnetic interactions (dipole–dipole and/or exchange ones) between grains on equilibrium^49–52 and magnetodynamic^53–55 response of a MMNP to an applied field. These works indicated fundamental similarities and profound differences between MMNPs and droplets in ferrofluid emulsions.^56–58 The similarities stem from the demagnetising effects, while the differences arise due to the spatial constraints. The latter is related to the fact that the MMNP grains are fixed in space and can only move together with the MMNP in contrast to ferroparticles in droplets that can freely diffuse inside. What remains unclear up to now is how the MMNPs self-assemble and what are the main parameters to alter, if one aims to control the clustering in these systems.

In this work, we will employ Langevin dynamics simulations and analytical theory in order to elucidate (i) MMNPs self-assembly in a wide range of magnetic interactions; (ii) the impact of MMNP self-assembly on the static magnetic susceptibility; in thermodynamic equilibrium.

Both questions raised above are particularly interesting as multicore magnetic particles appear to be similar to colloids with multiple patches,^59–63 but due to the possibility of the granular dipoles to reorient, multicores might exhibit a certain proximity to particles with mobile patches.^64–66 Moreover, recent experimental studies revealed a difference in self-assembly of magnetic nanoparticles and complex multicore colloids obtained by oil phase evaporation-induced self-assembly of hydrophobic nanoparticles.⁶⁷

The manuscript is organised as follows. Firstly, we analyse the initial susceptibility and find that, depending on the magnetic interactions between the grains, there are two distinct regimes: the system is dominated by demagnitising effects within individual MMNPs and remains very weakly susceptible; the system exhibits a steep increase in the initial susceptibility. Secondly, we show that this steep increase is caused by the MMNPs self-assembly. Within the clusters, MMNPs magnetise each other and become more correlated, albeit not as much as one could expect if the MMNPs were forming chains similar to their single-domain counterparts. Finally, we show that the reason for the grains to not form linear highly correlated structures is the formation of bridges between the MMNPs – clusters of correlated grains that go cross neighbouring MMNPs. Through bridges a MMNP might connect to more than two nearest neighbours. As a result, the clusters formed in the suspensions of MMNPs are more compact than those made by single-domain nanoparticles on the one hand, but, on the other, they are less compact than the bundle-like structures formed by magnetisable non-composite micron-sized magnetic particles.

2 Model of a multicore-based ferrofluid

We consider the system of N spherical MMNPs with the diameter D and volume v = πD³/6, dispersed in a three-dimensional nonmagnetic liquid medium with volume V.

Their volume fraction in the dispersion is

Each MMNP is filled with n_g spherical magnetic grains of diameter d (see Fig. 1). Positions of grains inside MMNPs are rigidly fixed. Grains are placed randomly and uniformly, without overlapping. Volume fraction of grains within each MMNP is

where v_g = πd³/6 is the grain volume. MMNPs in the medium are subject to both translational and rotational Brownian motion.

Fig. 1 Sketch of the investigated system.

The system is thermostated and has a constant temperature T, the external magnetic field is absent.

Grains are assumed to be single-domain. Each grain carries a magnetic moment [small mu, Greek, vector] of constant magnitude. We will make a simplifying assumption that the internal magnetic anisotropy energy of grains is negligible compared to thermal energy k_BT (k_B is the Boltzmann constant), and that magnetic moments can freely rotate under a combined effect of thermal fluctuations and a total dipolar field created by all other grains in the system. However, it is important to note that zero-field equilibrium properties of immobilized single-domain particles do not depend on the anisotropy energy, if the orientation distribution of their easy axes is random and uniform.^68–70 Such superparamagnetic systems sometimes referred to as “non-textured”.⁷¹ Therefore, the results obtained here for isotropic grains can also be extrapolated to non-textured MMNPs with anisotropic grains.

We account for the magnetic intergrain interactions by means of the conventional dipole–dipole pair potential

where [small mu, Greek, vector]_i and [small mu, Greek, vector]_j are the dipole moments of grains i and j, r⃑_ij = r⃑_i − r⃑_j is the displacement vector connecting their centers and r_ij = |r⃑_ij|, μ₀ is the magnetic permeability of vacuum. For characterizing magnetic interactions in our system we use the grain–grain dipolar coupling parameterthat is the characteristic energy scale of two adjacent grains with dipole moments, divided by k_BT.

In experiment, the number of grains per MMNP can vary from tens to hundreds.⁴⁸ In this work we will consider n_g = 20. This specific value was chosen as a practical compromise. On the one hand, it is big enough so that ensemble-averaged properties of the system can be analysed within a continuum-like approach.⁵¹ On the other, it is small enough so that we can simulate a suspension within a reasonable time frame. The grain packing fraction within MMNPs will be fixed to φ_g = 0.2. According to our packing fraction definition, eqn(2), it corresponds to a MMNP diameter of D ≃ 4.64d. For a typical value of d = 10 nm, D ≃ 46 nm. Thus, both the number of grains and the nanoparticle linear size lie within experimentally realistic ranges.⁷² It is worth noting, that for such MMNP configuration the majority of grains are actually positioned close to the nanoparticle surface. The fraction of near-surface grains can roughly be estimated as v_surf/v = 1 − (1 − 2d/D)³ ≃ 0.8, where v_surf is the volume of an MMNP outermost spherical shell of thickness d. This fact further strengthens the analogy between our system and patchy colloids pointed out in Section 1. An additional discussion on the potential importance of n_g and φ_g will be given in Section 3.2.2.

Notice, that in our model the grains are basically close-packed inside the multicore colloid and are not separated by any significant steric layers, similarly to the experimental realisation in ref. 73. In general the presence of the layers between the grains can lead to a decrease of the integrain interactions and as a result affect strongly the magnetic response of such particles to an applied magnetic field, as reported in the experimental study.⁷⁴

The main variables of the work are dipolar coupling constant 0 ≤ λ ≤ 10 and MMNP concentration 0.02 ≤ Φ ≤ 0.2. Langevin dynamics simulations will be used to determine equilibrium properties of MMNP suspensions within these parameter ranges. All the technical details of the method are given in Appendix A.

3 Results and discussion

3.1 Initial magnetic susceptibility

The key quantity that characterizes equilibrium zero-field properties of ferrofluids is the initial magnetic susceptibility χ (i.e., the initial slope of the magnetization curve). It can be calculated as⁷⁵Here, n_tot = n_g·N is the total number of dipoles (grains), M⃑ is the total magnetic moment of the system,〈⋯〉 denotes an ensemble average, 〈M⃑〉 is zero for our system in thermodynamic equilibrium due to symmetry consideration (there are no special directions in our system and all orientations are equally probable), χ_L is the so-called Langevin susceptibility,Langevin susceptibility characterizes the initial magnetostatic response of an ensemble of non-interacting dipoles (i.e., of an ideal superparamagnetic gas). Recently, we have suggested two simple analytical estimations for the susceptibility of MMNP ensembles.⁵² The first one takes into account dipole–dipole interactions between grains within individual MMNPs (including demagnetization effect), but neglects interactions between different MMNPs:Here, χ^g_L is the Langevin susceptibility, corresponding to a MMNP interior,eqn (8) was derived within the so-called modified mean-field theory.⁷⁶ The second equation adds a correction to take into account MMNPs' mutual magnetization (similarly to how it was done for ferroemulsions in ref. 56):

Simulation results for susceptibilities at different ferrofluid concentrations Φ are shown in Fig. 2. Values of χ are plotted against corresponding Langevin values χ_L (note the logarithmic scales). Using this representation it is easy to demonstrate the qualitative effect of dipole–dipole interactions on the system. If the susceptibility lies below the diagonal χ = χ_L, then the dipoles on average will try to prevent each other from aligning with an applied field, the overall magnetic response of the system will decrease. On the contrary, susceptibility values χ > χ_L indicate a reinforcing role of dipolar interactions. From this standpoint, simulation susceptibility curves in Fig. 2 demonstrate a non-trivial behaviour. In general, the Langevin ideal gas approximation only works in the limit χ_L ≪ 1 for all investigated samples. After that the susceptibility growth slows down and all the simulation curves fall below the diagonal. This behaviour is well captured by eqn (8) and (10). The predictions from two equations are close, the difference only becomes noticeable for the most concentrated samples. For Φ = 0.2, eqn (10) gives 25% larger susceptibilities than eqn (8), and it is actually much closer to simulation data up to χ_L ∼ 1. It indicates that at least for intermediate χ_L values the interactions between MMNPs are being taken into account correctly by eqn (10). However, both theories predict the susceptibility saturation at χ_L ≫ 1, which is mainly due to demagnetization effect within individual MMNPs.^51,52 Simulation curves, instead of reaching saturation, demonstrate a steep increase. This rapid change is happening at different χ_L for different samples. The inset in Fig. 2 shows the same susceptibility data plotted against dipolar coupling parameter λ. From these plots it is clear that the susceptibility increase is almost exclusively governed by an increase in λ. For our systems, it is always happening at λ ≳ 6, regardless the concentration value. The last important observation from Fig. 2 is that at λ ∼ 10 the susceptibility of concentrated samples with Φ ≥ 0.1 can actually reach values above Langevin ones. It means that the overall influence of magnetic interactions can switch from the suppression of the effective field–dipole coupling to its amplification.

Fig. 2 Initial magnetic susceptibility of a multicore-based ferrofluid vs. the corresponding Langevin susceptibility χ_L. Different colors correspond to different volume fractions of MMNPs (see colorbar). Circles show simulation results for n_g = 20 and φ_g = 0.2 (filled regions indicate 95% confidence intervals), dotted lines – theoretical predictions from eqn (8), Theory I, solid lines – theoretical predictions from eqn (10), Theory II. The dashed line corresponds to the susceptibility of an ideal superparamagnet with χ = χ_L. In the inset the same susceptibility values are plotted against the grain–grain dipolar coupling parameter λ.

We already know that at least in the presence of a strong applied field an anomalous increase in the magnetic response of MMNP ensembles is accompanied by a pronounced self-assembly.⁵² To understand whether the same is true in a zero-field case, a detailed cluster analysis was performed for one of the investigated concentration values, namely for a ferrofluid with Φ = 0.02. As the increase in susceptibility is mainly related to the change of λ, it is more important to vary the latter, rather than unnecessarily increase the set of sampling parameters by considering multiple concentrations and multiple magnetic coupling parameters at the same time.

3.2 Microstructure of a multicore-based ferrofluid

Looking at the upper row in Fig. 3, where the characteristic snapshots of the MMNP systems are presented, one can see how from an almost uniformly distributed in space (Fig. 3(a)) the dispersion exhibits a clear structural transformation if the value of λ increases and quite large clusters of MMNP are formed (Fig. 3(d)). Here, subfigures (a) and (e) are for λ = 6.25; (b) and (f) for λ = 7.5; (c) and (g) for λ = 8.75; and (d) and (h) for λ = 10. Even at the first glance one can notice that the topology of the clusters in Fig. 3 is notably different from those formed by dipolar hard or soft spheres. Compare, for example snapshots from works^77,78 to those in Fig. 3. Single-core magnetic particles mainly form structures with long linear segments and only few junctions. MMNPs, instead, seem to form more compact aggregates. If we zoom in the clusters, see the lower row of Fig. 3(e)–(h), some alignment of the dipoles in the granules of adjacent MNNPs can be observed. Thus, it is clear that the clustering in such system is occurring on two different spatial resolutions: MNNPs and individual granules. Below, as a consequence, we consider two different length-scales in order to explain equilibrium self-assembly of MMNPs.

Fig. 3 Representative simulation snapshots of the investigated systems for various values of dipole–dipole interactions: (a) and (e) λ = 6.25; (b) and (f) λ = 7.5; (c) and (g) λ = 8.75, (d) and (h) λ = 10. The bottom row is zoom in the clusters found in the dispersions from the upper row.

First, assuming each MMNP centre of mass to be a vertex of a graph, we identify the graph of connected components, meaning the clusters of MMNPs. The choice of a connected component is based on the definition of the edge connecting the two vertices in the graph. We consider two MMNPs as connected by an edge if the distance between their centres does not exceed 2^1/6D, and the total magnetic interaction between the two MMNPs is attractive, U_dd ≤ 0.

Second, we increase the resolution and look at the arrangement of individual grains in the previously defined clusters. We analyse how the grains of two adjacent MMNPs orient their magnetic dipoles in order to form the bond.

3.2.1 Clusters of MMNPs. In Fig. 4 we plot the cluster-size distributions. The vertical axis shows the probability of finding a cluster of a size shown along the horizontal axis. Each colour corresponds to a certain value of the magnetic coupling, as shown in the legend. This probability in case of single-core dipolar hard or soft spheres has an exponential form and should be linear if a log-scale is applied along the vertical axis as it is done in Fig. 4. However, as fits plotted with dotted lines show, only the decay λ = 6.25 might have been exponential (the third point should be taken with caution as there are very few clusters of the size three), all others are not. One can see that for the smallest value of λ, only about one per cent of dimers are there in the systems and the probability of finding a three particle cluster is negligible. For λ = 10, instead, only 60 per cent of MNNPs remain non-aggregated, while the rest of the system forms relatively large clusters.

Fig. 4 Cluster size distribution of a multicore-based ferrofluid. Values of dipole–dipole interactions are given in the legend. Dotted lines are guide for an eye to appreciate the deviations from the exponential decay observed for high values of λ.

As previously mentioned, not only the cluster size distributions, but also the topology of clusters formed by MMNPs seems to deviate from that of the dipolar hard or soft sphere particles. In order to verify this, for each MMNP, we calculate the number of nearest neighbours, i.e. the degree of each graph vertex. The histogram of MMNP degrees is shown in Fig. 5.

Fig. 5 Fraction of MMNPs having 0, 1, 2, 3, 4 or 5 neighbours. Values of dipole–dipole interactions are given in the legend.

Here, the same colour code as in Fig. 4 is used, although the results for λ = 6.25 are not provided as the majority of MMNPs are non-clustered. For λ = 7.5, only one per cent of MMNPs has two and more neighbours, meaning that only dimers form. For λ = 10, instead, more than 35 per cent of MMNPs have two and more neighbours, 30 per cent are having no neighbours. The remaining 35 per cent have one neighbour. The latter value shows that the number of open ends (vertices with degree one) typical for chain-like structures, or linear segments of branched clusters increases with magnetic interaction energy.

Even clearer picture of the cluster distribution and their topology can be obtained looking at Fig. 6. Here, typical graph representations of two simulation snapshots in equilibrium are presented. If λ = 6.25, Fig. 6(a), we observe one trimer, and seven dimers, all other MMNPs are nonaggregated. The system looks very different for λ = 10, Fig. 6(b). Here, we find several big connected components with vertices in them with degree two and higher. One should not be confused by stretched systems, as the representation of the graphs is simply optimised to occupy as few space as possible. In reality, those large clusters are 3D and rather compact, as the snapshot in Fig. 6(c) indicates. Here, for λ = 10, we show representative 3D structures of large clusters together with their graph representations, in order to avoid being confused by seemingly linear graphs.

Fig. 6 Typical graphs of the investigated MMNP dispersions for various values of dipole–dipole interactions: (a) λ = 6.25; (b) λ = 10. (c) Correspondence between graph representations and real clusters for λ = 10.

Previously, it was shown that the chain formation in single-core systems leads to the increase of the system magnetic response,⁷⁹ while the formation of rings reduces the susceptibility.⁸⁰ Below, we investigate how the new cluster topologies manifest themselves in the MMNP suspension magnetic properties.

3.2.2 Contribution of clusters into magnetostatic response. Combining the knowledge obtained form the inset of Fig. 2 and 4 it is obvious that the clustering of MMNPs and a susceptibility steep increase are happening at λ ≳ 6, so they must be correlated. To quantify the connection between two phenomena the following procedure was employed. Using the cluster analysis results, we were able to divide all MMNPs into two subgroups at each simulation step. The first subgroup contains non-clustered (single) MMNPs, let us denote their number as N_s. The second subgroup contains all MMNPs with non-zero degree and their number is denoted by N_c. Considering that the suspension is in a dynamic equilibrium, the values of N_c and N_s will fluctuate in time, but the sum of two numbers is fixed (N_s + N_c = N). We denote the total magnetic moment of all single MMNPs as M⃑_s, and the total magnetic moment of all clustered MMNPs as M⃑_c. Evidently, the total magnetisation of the suspension is the sum of the two, M⃑ = M⃑_s + M⃑_c. Using these definitions, we can rewrite the initial susceptibility eqn (5) as

χ = χ_s + χ_c + Δχ. (11)

Here, the expressions for the subgroups are:The last term in eqn (11) takes into account correlations between single particles and clusters:Of course, from the very nature of single particles, one can reasonably expect Δχ to be small.

Fig. 7(a) shows how χ_s and χ_c change with growing λ for Φ = 0.02. For reference, in the same figure we also provide the total susceptibility χ – both numerical and analytical – repeating the data from the inset of Fig. 2. Values of the correlation term are of the order of |Δχ| ∼ 10⁻³–10⁻² and are not shown. It is seen that up to the point where eqn (10) is still working, the susceptibility is mostly determined by single MMNPs, while the contribution of clusters is negligible (as well as their fraction in the system). After that, as the dipolar coupling parameter λ increases from 6.25 to 10, χ_s drops approximately by the factor of two and a half, that agrees well with Fig. 4, where the actual number of single particles drops from almost 200 at λ = 6.25 to approximately 80 at λ = 10.

Fig. 7 Zero-field magnetic characteristics of a multicore-based ferrofluid versus grain–grain dipolar coupling parameter λ. Φ = 0.02. Circles are simulation results averaged over the whole system, triangles correspond only to non-clustered (soliraty) MMNPs and squares correspond to clustered MMNPs. (a) Initial magnetic susceptibility, solid line is from eqn (10); (b) average normalized magnetic moment of a MMNP; (c) effective parameter of dipolar coupling between a pair of MMNPs.

In contrast to the susceptibility of single particles, the value of χ_c grows rapidly and at λ = 10 is twice as high as the theoretical saturation value χ ≃ 0.061 (as predicted by eqn (10)).

Another quantity that, along with susceptibility, can provide insights into magnetic properties of our system is the average magnetic moment of a MMNP. We calculate it as

where is the total instantaneous magnetic moment of the ith MMNP,Vectors [small mu, Greek, vector]⁽ⁱ⁾_j are the magnetic moments of the grains in the ith MMNP. The average magnetic moment m is normalized, its maximal value m = 1 corresponds to the magnetic saturation. If interactions between magnetic moments of grains are weak, m must be close to n_g^−1/2 ≃ 0.22. We can also calculate average magnetic moments separately for two particle subgroups:where and denote summation over single and clustered MMNPs, respectively.

In Fig. 7(b) we plot dependencies of magnetic moments m, m_s and m_c on λ. Considering that the MMNPs form a statistically relevant number of clusters only for λ ≥ 5, at lower coupling constants the number of clustered particles N_c is too small to reliably estimate m_c from eqn (18). Initially, magnetic moments are decreasing with λ. This is the result of the demagnetisation effect – for a localised collection of interacting dipoles surrounded by a non-magnetic medium it is always energetically more favorable to minimise its total magnetic moment.⁸¹ However, while for single particles this decrease continues up to λ = 10, the dependence of m is found to be non-monotonic. At λ ≳ 6, as the fraction and the average size of clusters start to grow, the magnetic moment of clustered particles turns into an increasing function as well. MMNPs within clusters, therefore, help each other to overcome the demagnetisation effect and become more susceptible to the applied field.

Finally, knowing numerical values for MMNPs magnetic moments, we can draw a comparison with dipolar spheres, which are used for modelling of traditional single-core-based ferrofluids. For each value of λ, instead of an ensemble of MMNPs, let us consider an equivalent ensemble of dipolar spheres of diameter D with a point-like dipole in its center. The magnitude of this magnetic moment is constant and equals to where m is calculated from eqn (15). Then similarly to eqn (4) we can introduce an effective dipolar coupling parameter for a pair of MMNPs as

A notable thing in eqn (19) is that it suggests that interactions between MMNPs must strongly depend on both number of grains n_g and their packing fraction φ_g. We expect that changing these parameters can change a critical value of λ, at which the self-assembly starts to influence system magnetic response (in our case, this is λ ≃ 6). A detailed inquiry into the importance of n_g and φ_g is left for future studies.

Values of Λ for n_g = 20 and φ_g = 0.02 are shown in Fig. 7(c). We can see that it actually remains quite low, Λ < 1. If we estimate the coupling parameter independently for two particle subgroups (Λ_s = λn_gφ_gm_s² and Λ_c = λn_gφ_gm_c²), it will be somewhat larger for clustered particles and can reach Λ_c = 1.2. But we know that for real dipolar spheres a pronounced self-assembly in the absence of a field only takes place at Λ ≥ 4^77,82! Therefore, the crude analogy with dipolar spheres is insufficient to explain the self-assembly mechanism in multicore-based ferrofluids. An obvious difference between two systems is that the total magnetic moment of a multicore particle is not point-like, it is distributed across the particle volume. Presumably, correlations between individual grains within interacting MMNPs must be explicitly taken into account in order to explain the observed behavior. In fact, below we reveal the mechanism that keeps MMNPs together even though the effective Λ between them is very low.

3.2.3 Bridges between multicore nanoparticles. Increasing the resolution of the analysis, in this subsection, we study the internal structure of clusters, namely, the correlations between the grains inside MMNPs. Looking closer at Fig. 3(e)–(h), one can notice that individual grains from different MMNPs cooperate and form their own internal clusters. Such type of clusters, consisting of grains belonging to different MMNPs, we call bridges and provide a more detailed snapshot in Fig. 8.

Fig. 8 Snapshot of a MMNP cluster. Grains on the borders between different MMNPs, forming bridges, are encircled.

In order to analyze the structure and number of bridges, we once again carried out a cluster analysis. Firstly, we localize all MMNP clusters as connected components (see, explanations above) and sort them into groups according to the number of MMNPs in them. Secondly, for each such element of each subgroup, we construct a graph whose vertices are the individual grains and find connected components within these new structures using two conditions (i) two grains are connected if the distance between them is smaller than 2^1/6d and the dipolar interaction between them is negative; (ii) the connected component is a cluster of interest if it contains vertices from different MMNPs. Finally, any property of bridges that is calculated is averaged over the MMNP clusters of a given size.

In Fig. 9, we present the probability of finding a certain number of bridges within the cluster of two (a), three (b) and four (c) MMNPs. If λ = 7.5 the only reliable statistics can be obtained for dimers. Here, half of the clusters have at least one bridge, approximately 12 per cent have two bridges holding the particles together. We believe, that 35 per cent of clusters without bridges is related to a very restrictive distance criteria that we apply to detect clusters. In fact, with growing dipolar coupling one sees that the number of dimers without connecting bridges vanishes, while the probability of finding two bridges is increasing. Three particle clusters are normally connected with two-three bridges. Four-particle clusters are forming frequently enough only for λ = 10. Here within the 10 per cent margin one can find up to five bridges.

Fig. 9 Number of bridges in MMNP clusters of different sizes: (a) cluster size equals 2, (b) cluster size equals 3, (c) cluster size equals 4. Values of dipole–dipole interactions are given in the legend.

In Fig. 10 one can see that the bridges are rather long. With only slightly decreasing probability, one can find bridges with length from two to seven particles in dimers Fig. 10(a). Interestingly enough, comparing Fig. 10(b) and (c), where the bridge length distributions are plotted for three- and four-particle clusters respectively, to the histogram for dimers in (a), one concludes that the difference is barely to be found. In other words, the length of the bridges is cluster-size independent and it is only weakly affected by λ. This suggests that the bridges that we find reach their limiting length determined by spatial constraints of the MMNP size. In fact, if grains follow MMNP diameter, only 4 of them can fit inside. Taking this into account, it is particularly curious to elucidate what is a typical bridge topology and if it depends on the bridge length.

Fig. 10 Bridge size in MMNP clusters of different values: (a) cluster size equals 2, (b) cluster size equals 3, (c) cluster size equals 4. Values of dipole–dipole interactions are given in the legend.

The first step to define the topology of the bridges is to calculate the degree of the grains in bridges. The results are collected in Fig. 11. Had the bridges been chains, there would have been only grains with degrees one and two. If, for example, the chain-like bridge was built of six grains, roughly 33 per cent would have degree one and 67 per cent degree two. The longer the chain, the lower the fraction of free ends in comparison to the grains with degree two. Here, however independently from the values of λ and the cluster size, 35–40 per cent of grains have degree one, 40–45 per cent degree two, around 15 per cent of grains in bridges have degree of three and five per cent of grains have four different neighbours.

Fig. 11 Fraction of grains having a given number of neighbours in bridges for various dipole–dipole interactions: (a) λ = 7.5, (b) λ = 8.75, (c) λ = 10. The averaging is performed among MMNP clusters of different sizes as shown in the legend.

Fig. 12(a)–(c) show the total magnetic moment of the bridge, normalized by the magnitude of the grain magnetic moment and bridge length. It can be seen that for all values of λ and cluster sizes, there is a bimodal distribution of the magnetic moment: the first peak in the region of 0.2–0.3 and there is a main maximum near unity. The latter corresponds to the linear structure of the bridge, whereas the peak at 0.2 indicates some bend, possibly branched or closed, configuration. Importantly, the bend structures are not ideal rings, as the magnetisation is not that low as one would expect for a mildly perturbed ring. Fig. 12(d) shows the dependence of the average bridge magnetic moment on the bridge length. It confirms the presence of long banded or nearly closed structures: short bridges have a very high magnetic moment, while the second histogram maximum in Fig. 12(a)–(c), shown by a vertical line, corresponds to the bridges longer than 10. It means that grains with degrees three and four are the bending points. Long banded bridges often start and end in the same MMNP, while the part in the fold belongs to the neighbouring MMNP in the cluster. Only few very long non-banded bridges were found in simulations that would percolate through three or more grains.

Fig. 12 Average magnetic moment of bridges in MMNP clusters of different sizes: (a) cluster size equals 2, (b) cluster size equals 3, (c) cluster size equals 4. (d) The magnetic moment of bridges depending on the bridge size. Values of dipole–dipole interactions are given in the legend.

The formation of various bridges allows to increase the number of nearest neighbours a MMNP can have, making it similar to colloids with mobile multiple patches.⁸³

4 Conclusions

In this work, using Langevin dynamics computer simulations accompanied by analytical calculations, we investigated diluted suspensions of multicore magnetic nanosized particles (MMNPs) composed by grains, whose positions are fixed within the particle body, but their magnetic moments are free to rotate, corresponding to the grains with negligibly low magnetic anisotropy.

Self-assembly of MMNPs and its impact on the zero-field magnetostatic response of the MMNPs suspensions were addressed at two different scales.

Firstly, it was shown that MMNPs started forming clusters at the values of saturation magnetisation much higher than their point-dipolar counterparts. Instead of forming inherent to dipolar spheres chains or rings whose size distribution decays exponentially with the cluster size, for higher interaction strength, the clusters formed by MMNPs are more compact, and their sizes are not distributed exponentially. In fact, the total average magnetic moment of a particle is rather small, and for the conditions in which self-assembly does not take place, the static susceptibility of MMNP suspensions is found to be below Langevin one and to increase only slightly with growing magnetic coupling parameter. At the value of λ at which the self-assembly ushers in, the susceptibility of the system qualitatively changes, starting growing very fast and going above the susceptibility of an ideal superparamagnetic gas predicted by Langevin theory.

Secondly, in order to understand the structure of the clusters formed by MMNPs, we looked at the interactions of the individual grains in neighbouring MMNPs. We discovered that clustering occurs by means of building grain bridges: grains of kissing MMNPs reorient so that they can form clusters. The length of those clusters (bridges) is found to be only weakly dependent on the magnetic interaction strength if the latter is high enough. One can find bridges whose length varies between 3–4 grains to 7–11 of them. Short bridges are mainly linear and have a relatively large magnetic moment. Long bridges instead have banded, u-like shape that is nearly closed with rather small total magnetic moment compensated by their folded structure.

Bridge formations is responsible for the multiple bonds that MMNPs can form, attaching more than two nearest neighbours.

This study reveals a mechanism of MMNPs self-assembly that makes this system clearly different from both single-domain nanoparticles and magnetisable micron-sized colloids. As the next step, we would address the dynamics of MMNPs.

Even though, the direct discussion of the dynamical properties of multicore suspensions cannot be based on the results presented in this paper, we cannot but mention the potential impact of structural transitions on the dynamic susceptibility of these systems. We assume that the formation of bridges as a mean for multicore particles to self-assemble will lead to a significant broadening of the dynamic spectra. At the moment we are working on the generalisation of our approach to capture off-equilibrium dynamics.

Conflicts of interest

There are no conflicts to declare.

Appendices

A simulation details

All the results reported here are from molecular dynamics simulations performed in the software package ESPResSo 4.1.4.⁸⁴ The “virtual sites” feature of ESPResSo is used to take into account the granular nature of MMNPs.⁸⁵ Each magnetic grain is treated as a virtual site, connected to the center of the corresponding “host” MMNP. The position of each virtual site is rigidly fixed within the body reference frame of the host. All forces and torques acting on virtual sites are propagated back to the host on every simulation time step. Langevin thermostat technique is used to achieve thermodynamic equilibrium. Formally, Langevin equations for translation and rotational motion of a given MMNP can be written aswhere asterisk denotes reduced quantities, d is used as a unit of length, grain mass p is used as a unit of mass and k_BT is used as a unit of energy. Thus, and are the reduced linear and angular velocities of a MMNP, respectively, P* = P/p is the reduced MMNP mass, J* = J/pd² is its reduced moment of inertia, is the total reduced force acting on a given MMNP, f⃑^C is the central force acting on a MMNP due to steric interactions, f⃑^g_j is the total force acting on the jth virtual site within MMNP, , r⃑^C_j is the virtual site position relative to the MMNP center of mass, [small tau, Greek, vector]^g_j is the total torque on the jth virtual site, and are the reduced translational and rotational friction coefficients, and are the random thermal force and torque on a MMNP, respectively. They have zero mean values and are connected to and via standard fluctuation-dissipation relationships.⁸⁶ The reduced time is .

While the “virtual sites” feature is used to calculate the positions of superparamagnetic grains, the equilibrium orientations of their magnetic moments are determined independently from Langevin equations identical in form to eqn (21):

Here, is the reduced magnetic torque acting on a given grain, H⃑_dd is the sum all dipolar fields in the grain centre, dimensionless quantities , , and have the same meaning as their counterparts from eqn (21).

To simulate the steric repulsion between MMNPs, the Weeks–Chandler–Andersen (WCA) pair potential is used:⁸⁷

where R is the distance between centers of two MMNPs. Additionally, a similar potential is imposed between each pair of grains in the system:Introduction of this potential was required to avoid numerical instabilities arising in rare cases when the dipole–dipole interaction between two grains from two different MMNPs was large enough to overcome a soft core repulsion modelled by eqn (23) (at λ ∼ 10). This approach allowed us to use a larger time step.

All simulations are performed for N = 200 MMNPs placed inside a cubic box with a side length of . 3D periodic boundary conditions are imposed on a box. The torques due to long-range dipole–dipole interactions are computed using the dipolar P³M algorithm.⁸⁸ Other simulation parameters are P* = N_g = 20, J* = 0.1P*(D/d)² ≃ 43, . The simulation time step is Δt* = 0.002. Typically, the first 5 × 10⁵ time steps are used for the system equilibration, and the subsequent production run lasts for at least another 10 × 10⁵ time steps. For every combination of Φ and λ at least five independent runs are performed.

Albeit high efficiency, the choice of such dimensionless parameters imposes a certain “limitation”: the dynamics of our system is not realistic, and we cannot extract any off-equilibrium properties of the system. Although very interesting, the life-time of the clusters and dynamic response of the multicores lay beyond the scope of this study and require more frequent sampling with a smaller time-step, not compatible with reaching equilibrium in a feasible time.

Acknowledgements

This research has been supported by the Russian Science Foundation Grant No. 19-72-10033. S. S. K. was partially supported by Project SAM P 33748. Computer simulations were performed at the Ural Federal University cluster.

Notes and references

1. R. Tietze, J. Zaloga, H. Unterweger, S. Lyer, R. P. Friedrich, C. Janko, M. Pöttler, S. Dürr and C. Alexiou, Biochem. Biophys. Res. Commun., 2015, 468, 463–470.
2. C. Ménager, O. Sandre, J. Mangili and V. Cabuil, Polymer, 2004, 45, 2475–2481.
3. S. Backes, M. U. Witt, E. Roeben, L. Kuhrts, S. Aleed, A. M. Schmidt and R. von Klitzing, J. Phys. Chem. B, 2015, 119, 12129–12137.
4. P. Mandal, S. Maji, S. Panja, O. P. Bajpai, T. K. Maiti and S. Chattopadhyay, New J. Chem., 2019, 43, 3026–3037.
5. B. Sung, M.-H. Kim and L. Abelmann, Bioeng. Trans. Med., 2021, 6, e10190.
6. Y. Cao, Z. Mao, Y. He, Y. Kuang, M. Liu, Y. Zhou, Y. Zhang and R. Pei, ACS Appl. Mater. Interfaces, 2020, 12, 26973–26981.
7. F. Gao, X. Wu, D. Wu, J. Yu, J. Yao, Q. Qi, Z. Cao, Q. Cui and Y. Mi, Colloids Surf., A, 2020, 587, 124363.
8. C. Biglione, J. Bergueiro, S. Wedepohl, B. Klemke, M. C. Strumia and M. Calderón, Nanoscale, 2020, 12, 21635–21646.
9. I. Hilger and W. A. Kaiser, Nanomedicine, 2012, 7, 1443–1459.
10. S. Dürr, C. Janko, S. Lyer, P. Tripal, M. Schwarz, J. Zaloga, R. Tietze and C. Alexiou, Nanotechnol. Rev., 2013, 2, 395–409.
11. E. A. Périgo, G. Hemery, O. Sandre, D. Ortega, E. Garaio, F. Plazaola and F. J. Teran, Appl. Phys. Rev., 2015, 2, 041302.
12. E. C. Abenojar, S. Wickramasinghe, J. Bas-Concepcion and A. C. S. Samia, Prog. Natl. Sci.: Mater. Int., 2016, 26, 440–448.
13. E. Resler Jr and R. Rosensweig, AIAA J., 1964, 2, 1418–1422.
14. E. M. Furst, C. Suzuki, M. Fermigier and A. P. Gast, Langmuir, 1998, 14, 7334–7336.
15. C. Goubault, P. Jop, M. Fermigier, J. Baudry, E. Bertrand and J. Bibette, Phys. Rev. Lett., 2003, 91, 260802.
16. J. J. Benkoski, S. E. Bowles, R. L. Jones, J. F. Douglas, J. Pyun and A. Karim, J. Polym. Sci., Part B: Polym. Phys., 2008, 46, 2267–2277.
17. Z. Zhou, G. Liu and D. Han, ACS Nano, 2009, 3, 165–172.
18. J. J. Benkoski, J. L. Breidenich, O. M. Uy, A. T. Hayes, R. M. Deacon, H. B. Land, J. M. Spicer, P. Y. Keng and J. Pyun, J. Mater. Chem., 2011, 21, 7314–7325.
19. A. A. Kuznetsov, J. Magn. Magn. Mater., 2019, 470, 28–32.
20. R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone and J. Bibette, Nature, 2005, 437, 862–865.
21. É. Bereczk-Tompa, F. Vonderviszt, B. Horváth, I. Szalai and M. Pósfai, Nanoscale, 2017, 9, 15062–15069.
22. S. Kralj and D. Makovec, ACS Nano, 2015, 9, 9700–9707.
23. Y. Xiong, Z. Lin, D. Mostarac, B. Minevich, Q. Peng, G. Zhu, P. A. Sánchez, S. Kantorovich, Y. Ke and O. Gang, Nano Lett., 2021, 21, 10547–10554.
24. D. Mostarac, P. A. Sánchez and S. Kantorovich, Nanoscale, 2020, 12, 13933–13947.
25. J. Witte, T. Kyrey, J. Lutzki, A. M. Dahl, J. Houston, A. Radulescu, V. Pipich, L. Stingaciu, M. Kühnhammer, M. U. Witt, R. von Klitzing, O. Holderer and S. Wellert, Soft Matter, 2019, 15, 1053–1064.
26. I. S. Novikau, P. A. Sánchez and S. S. Kantorovich, J. Mol. Liq., 2020, 307, 112902.
27. D. Eberbeck, C. L. Dennis, N. F. Huls, K. L. Krycka, C. Gruttner and F. Westphal, IEEE Trans. Magn., 2013, 49, 269–274.
28. F. Ludwig, O. Kazakova, L. F. Barqun, A. Fornara, L. Trahms, U. Steinhoff, P. Svedlindh, E. Wetterskog, Q. A. Pankhurst and P. Southern, et al. , IEEE Trans. Magn., 2014, 50, 1–4.
29. P. Bender, D. Honecker and L. Fernández Barqun, Appl. Phys. Lett., 2019, 115, 132406.
30. S. B. Trisnanto and Y. Takemura, J. Appl. Phys., 2021, 130, 064302.
31. V. Socoliuc, M. Avdeev, V. Kuncser, R. Turcu, E. Tombácz and L. Vekas, Nanoscale, 2022, 14, 4786–4886.
32. A. Bakandritsos, A. Papagiannopoulos, E. N. Anagnostou, K. Avgoustakis, R. Zboril, S. Pispas, J. Tucek, V. Ryukhtin, N. Bouropoulos, A. Kolokithas-Ntoukas, T. A. Steriotis, U. Keiderling and F. Winnefeld, Small, 2012, 8, 2381–2393.
33. S. Lyer, R. Tietze, H. Unterweger, J. Zaloga, R. Singh, J. Matuszak, M. Poettler, R. P. Friedrich, S. Duerr, I. Cicha, C. Janko and C. Alexiou, Nanomedicine, 2015, 10, 3287–3304.
34. S. Ziemian, N. Löwa, O. Kosch, D. Bajj, F. Wiekhorst and G. Schütz, Nanomaterials, 2018, 8, 180.
35. S. Dutz, M. Kettering, I. Hilger, R. Müller and M. Zeisberger, Nanotechnology, 2011, 22, 265102.
36. C. Blanco-Andujar, D. Ortega, P. Southern, Q. A. Pankhurst and N. T. K. Thanh, Nanoscale, 2015, 7, 1768–1775.
37. G. Hemery, C. Genevois, F. Couillaud, S. Lacomme, E. Gontier, E. Ibarboure, S. Lecommandoux, E. Garanger and O. Sandre, Mol. Syst. Des. Eng., 2017, 2, 629–639.
38. P. Bender, J. Fock, M. Hansen, L. Bogart, P. Southern, F. Ludwig, F. Wiekhorst, W. Szczerba, L. Zeng and D. Heinke, et al. , Nanotechnology, 2018, 29, 425705.
39. J. Ko and H. B. Lim, Anal. Chem., 2014, 86, 4140–4144.
40. A. Jeong and H. Lim, Talanta, 2018, 178, 916–921.
41. J. Nowak, F. Wiekhorst, L. Trahms and S. Odenbach, J. Phys.: Condens. Matter, 2014, 26, 176004.
42. D. Zablotsky, S. Kralj, G. Kitenbergs and M. M. Maiorov, J. Non-Newtonian Fluid Mech., 2020, 278, 104248.
43. A. O. Ivanov and A. Zubarev, Materials, 2020, 13, 3956.
44. D. Borin, A. Zubarev, D. Chirikov, R. Müller and S. Odenbach, J. Magn. Magn. Mater., 2011, 323, 1273–1277.
45. D. Y. Borin, A. Y. Zubarev, D. N. Chirikov and S. Odenbach, J. Phys.: Condens. Matter, 2014, 26, 406002.
46. D. Zablotsky, S. Kralj and M. M. Maiorov, Colloids Surf., A, 2020, 603, 125079.
47. J. Wells, O. Kazakova, O. Posth, U. Steinhoff, S. Petronis, L. K. Bogart, P. Southern, Q. Pankhurst and C. Johansson, J. Phys. D, 2017, 50, 383003.
48. L. Gutiérrez, R. Costo, C. Grüttner, F. Westphal, N. Gehrke, D. Heinke, A. Fornara, Q. A. Pankhurst, C. Johansson, S. Veintemillas-Verdaguer and M. P. Morales, Dalton Trans., 2015, 44, 2943–2952.
49. V. Schaller, G. Wahnström, A. Sanz-Velasco, S. Gustafsson, E. Olsson, P. Enoksson and C. Johansson, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 092406.
50. V. Schaller, G. Wahnström, A. Sanz-Velasco, P. Enoksson and C. Johansson, J. Magn. Magn. Mater., 2009, 321, 1400–1403.
51. A. A. Kuznetsov, Phys. Rev. B, 2018, 98, 144418.
52. A. A. Kuznetsov, E. V. Novak, E. S. Pyanzina and S. S. Kantorovich, J. Mol. Liq., 2022, 119373.
53. A. Weddemann, A. Auge, D. Kappe, F. Wittbracht and A. Hütten, J. Magn. Magn. Mater., 2010, 322, 643–646.
54. N. A. Usov, O. N. Serebryakova and V. P. Tarasov, Nanoscale Res. Lett., 2017, 12, 489.
55. P. Ilg, Phys. Rev. B, 2017, 95, 214427.
56. A. O. Ivanov and O. B. Kuznetsova, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 041405.
57. I. M. Subbotin, Magnetohydrodynamics, 2018, 54, 131–136.
58. I. M. Subbotin, J. Magn. Magn. Mater., 2020, 502, 166524.
59. A. B. Pawar and I. Kretzschmar, Macromol. Rapid Commun., 2010, 31, 150–168.
60. H. Zhang and M. Widom, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 8951.
61. Z. Gong, T. Hueckel, G.-R. Yi and S. Sacanna, Nature, 2017, 550, 234–238.
62. F. Sciortino and E. Zaccarelli, Curr. Opin. Solid State Mater. Sci., 2011, 15, 246–253.
63. E. Bianchi, R. Blaak and C. N. Likos, Phys. Chem. Chem. Phys., 2011, 13, 6397–6410.
64. S. Angioletti-Uberti, P. Varilly, B. M. Mognetti and D. Frenkel, Phys. Rev. Lett., 2014, 113, 128303.
65. A. McMullen, M. Holmes-Cerfon, F. Sciortino, A. Y. Grosberg and J. Brujic, Phys. Rev. Lett., 2018, 121, 138002.
66. P. A. Sánchez, A. Caciagli, S. S. Kantorovich and E. Eiser, J. Mol. Liq., 2023, 382, 121895.
67. P. Qiu, C. Jensen, N. Charity, R. Towner and C. Mao, J. Am. Chem. Soc., 2010, 132, 17724–17732.
68. C. P. Bean and J. D. Livingston, J. Appl. Phys., 1959, 30, S120–S129.
69. E. A. Elfimova, A. O. Ivanov and P. J. Camp, Nanoscale, 2019, 11, 21834–21846.
70. A. O. Ivanov and F. Ludwig, Phys. Rev. E, 2020, 102, 032603.
71. F. G. Silva, R. Aquino, F. A. Tourinho, V. Stepanov, Y. L. Raikher, R. Perzynski and J. Depeyrot, J. Phys. D: Appl. Phys., 2013, 46, 285003.
72. S. Dutz, IEEE Trans. Magn., 2016, 52, 1–3.
73. A. Bunge, A. S. Porav, G. Borodi, T. Radu, A. Pîrnau, C. Berghian-Grosan and R. Turcu, J. Mater. Sci., 2019, 54, 2853–2875.
74. C. L. Dennis, K. L. Krycka, J. A. Borchers, R. D. Desautels, J. Van Lierop, N. F. Huls, A. J. Jackson, C. Gruettner and R. Ivkov, Adv. Funct. Mater., 2015, 25, 4300–4311.
75. S. W. De Leeuw, J. W. Perram and E. R. Smith, Annu. Rev. Phys. Chem., 1986, 37, 245–270.
76. A. O. Ivanov and O. B. Kuznetsova, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 041405.
77. L. Rovigatti, J. Russo and F. Sciortino, Soft Matter, 2012, 8, 6310–6319.
78. L. Rovigatti, J. Russo and F. Sciortino, Phys. Rev. Lett., 2011, 107, 237801.
79. V. Mendelev and A. Ivanov, J. Magn. Magn. Mater., 2005, 289, 211–214.
80. S. Kantorovich, A. O. Ivanov, L. Rovigatti, J. M. Tavares and F. Sciortino, Phys. Rev. Lett., 2013, 110, 148306.
81. C. Kittel, Introduction to solid state physics, Wiley, 2004.
82. A. O. Ivanov, Z. Wang and C. Holm, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 031206.
83. S. Angioletti-Uberti, P. Varilly, B. M. Mognetti and D. Frenkel, Phys. Rev. Lett., 2014, 113, 128303.
84. F. Weik, R. Weeber, K. Szuttor, K. Breitsprecher, J. de Graaf, M. Kuron, J. Landsgesell, H. Menke, D. Sean and C. Holm, Eur. Phys. J.: Spec. Top., 2019, 227, 1789–1816.
85. A. Arnold, O. Lenz, S. Kesselheim, R. Weeber, F. Fahrenberger, D. Roehm, P. Košovan and C. Holm, Meshfree methods for partial differential equations VI, Springer, 2013, pp. 1–23.
86. W. T. Coffey, Y. P. Kalmykov and J. T. Waldron, The Langevin equation: with applications to stochastic problems in physics, chemistry and electrical engineering, World Scientific, Singapore, 2004.
87. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
88. J. J. Cerdà, V. Ballenegger, O. Lenz and C. Holm, J. Chem. Phys., 2008, 129, 234104.

---