Universal behavior of dense clusters of magnetic nanoparticles

A detailed numerical simulation of quasistatic hysteresis loops of dense clusters of interacting magnetic nanoparticles is carried out. Both clusters of magnetically soft and magnetically hard nanoparticles are considered. The clusters are characterized by an average particle diameter D, the cluster radius Rc, the particle saturation magnetization Ms, and the uniaxial anisotropy constant K. The number of particles in the cluster varies between Np = 30 - 120. The particle centers are randomly distributed within the cluster, their easy anisotropy axes being randomly oriented. It is shown that a rare assembly of identical random clusters of magnetic nanoparticles can be characterized by two dimensionless parameters: 1) the relative strength of magneto-dipole interaction, K/Ms^2, and the average particle concentration within the cluster, {\eta} = VNp/Vc. Here V is the nanoparticle volume, and Vc is the volume of the cluster, respectively. In the strong interaction limit, Ms*{\eta}/Ha>>1, where Ha = 2K/Ms is the anisotropy field, the ultimate hysteresis loops of dilute assemblies of clusters have been constructed. In the variables (M/Ms, H/Ms) these hysteresis loops depend only on the particle volume fraction {\eta}. In the weak interaction limit, Ms*{\eta}/Ha<<1, the assembly hysteresis loops in the variables (M/Ms, H/Ha) are close to the standard Stoner-Wohlfarth hysteresis loop.


I. INTRODUCTION
Dense assemblies of magnetic nanoparticles show very rich and complex behavior because their properties depend on several important factors, such as the geometrical structure of the assembly, the nature of magnetic anisotropy of isolated nanoparticles, the presence of exchange and magneto-dipole interaction among the nanoparticles, etc. The long -range character of the inter-particle magneto-dipole interaction is the main problem that complicates significantly theoretical investigation of the dense nanoparticle assemblies. For comparison, the properties of assemblies of non interacting single-domain nanoparticles have been studied [1][2][3][4][5][6][7][8][9][10] in detail for nanoparticles with uniaxial, cubic and combined types of magnetic anisotropy. In addition, the behavior of dilute assemblies of superparamagnetic nanoparticles has been investigated depending on the ambient medium temperature and the value of the phenomenological damping constant [11 -15]. However, theoretical study of dense assemblies of magnetic nanoparticles is more difficult as for an assembly of N p nanoparticles N p 2 pair interactions must be taken into account. Therefore, the computational complexity increases rapidly with the number of the particles in the assembly.
In spite of various interesting results obtained the understanding of the properties of dense assemblies of magnetic nanoparticles is still incomplete, since in general it is necessary to investigate the behavior of an assembly of interacting nanoparticles in solid matrix and viscous liquid, at moderate and high temperatures, under the influence of the quasi-static and alternating external magnetic field, etc. In this paper the quasi-static hysteresis loops of an assembly of magnetic nanoparticles distributed in a solid matrix are studied using numerical simulation under the following assumptions: 1) spherical nanoparticles of a radius R are characterized by the saturation magnetization M s and uniaxial magnetic anisotropy constant K; 2) the temperature of the medium is much lower than the particle blocking temperature T b , so that the thermal fluctuations of the particle magnetic moments can be neglected; 3) the particles of the assembly are contained in the quasi-spherical clusters of the given radius R cl ; 4) each cluster contains approximately the same number of nanoparticles N p >> 1. The assembly of clusters is sufficiently sparse, so that the magnetostatic interaction of nanoparticles of different clusters can be neglected. It is also assumed that the easy anisotropy axes of nanoparticles are oriented randomly, and the centers of nanoparticles occupy random positions within the cluster volume. However, the nanoparticles do not touch each other, so that the exchange interaction between neighboring nanoparticles is absent.
In spite of the constraints assumed, this model seems to be able to describe adequately the behavior the experimentally investigated assemblies of clusters with strong magneto-dipole interaction [35][36][37][38][39][40]. Note that even with the simplifications adopted, there are 5 independent parameters of the model, i.e. the saturation magnetization M s , the uniaxial anisotropy constant K, the particle radius R, the cluster radius R cl , and the total number of the nanoparticles in the cluster N p . In this paper, based on the detailed numerical simulations, a unified description of the quasi-static hysteresis loops of the assembly of dense clusters of nanoparticles is obtained. The hysteresis loops of the assembly are averaged over a sufficient number (~ 20 -30) of random cluster realizations. The total number of particles in the whole assembly is about of several thousands. This is sufficient to obtain statistically convincing results.

II. Basic equations
Dynamics of the unit magnetization vector α i of the i-th nanoparticle of the cluster is described by the LLG equation where κ is the phenomenological damping constant and γ is the gyromagnetic ratio. The effective magnetic field acting on a separate nanoparticle can be calculated as a derivative of the total cluster energy ( Here 6 3 D V π = is the nanoparticle volume, D = 2R is the particle diameter. The total magnetic energy of the cluster W = W a + W Z + W m is a sum of the magnetocrystalline anisotropy energy W a , Zeeman energy W Z of the particles in the external uniform magnetic field H, and the energy of mutual magneto-dipole interaction of the particles W m .
For nanoparticles of nearly spherical shape with uniaxial type of magneto-crystalline anisotropy the magneto-crystalline anisotropy energy is given by where e i is the orientation of the easy anisotropy axis of i-th particle of the cluster. Zeeman energy W Z of the cluster in the external uniform magnetic field is given by Next, for spherical uniformly magnetized nanoparticles the magnetostatic energy of the cluster can be represented as the energy of the point interacting dipoles. Let r i be the coordinates of the particle centers within the cluster. Then the magneto-dipole interacting energy reads where n ij is the unit vector along the line connecting the centers of i-th and j-th particles, respectively. Thus, the effective magnetic field acting on the i-th nanoparticle of the cluster is given by where H a = 2K/M s is the particle anisotropy field.
We start the calculation of a quasistatic hysteresis loop of a cluster at a sufficiently high value of the external uniform magnetic field sufficient for assembly saturation. Therefore, in the initial state the unit magnetization vectors of all particles are set in the direction parallel to the applied magnetic field. At the next step the magnetic field decreases by a small value, dH = 1 -2 Oe, and evolution of the initial magnetization state is calculated according to Eqs. (1) - (5). For a cluster consisting of several tens of nanoparticles it is necessary to make usually 10 4 -10 5 iterations with a sufficiently small time increment to find stable magnetization state of a cluster at a given value of the external magnetic field. In accordance with the Eq. (1), the final magnetization state is assumed to be stable under the condition which means that the unit magnetization vector of each nanoparticle is parallel with a sufficiently high accuracy to the effective magnetic field acting on this nanoparticle. In this manner, for a given set of initial material and geometrical parameters, i.e. K, M s , R, R cl and N p , one can obtain a quasi-static hysteresis loop of a nearly spherical cluster of nanoparticles, which depend, generally speaking, on the set of the coordinates of the nanoparticle centers {r i }, as well as on the set the orientations {e i } of the particle easy anisotropy axes. However, the calculations show that in the limit N p >> 1 the hysteresis loops obtained for different realizations of random variables {r i } and {e i } differ only slightly from each other.
Moreover, if one assume an assembly of noninteracting identical clusters, and calculate a common hysteresis loop of the assembly averaged over a sufficiently large number of random cluster realizations, one obtains a non-random hysteresis loop. The latter characterizes the behaviour of the assembly of noninteracting random nanoparticle clusters. Calculations show that in the limit N p >> 1 the averaged hysteresis loop of cluster assembly has a rather small dispersion even being averaged over 20 -30 independent realizations of random clusters with the fixed values of the initial parameters K, M s , R, R cl and N p .
The question naturally arises, how to describe these non-random hysteresis loops and to establish an actual set of dimensionless parameters which determine the averaged hysteresis loop behaviour. Note that the average value of the effective magnetic field, equation (6), acting on a typical particle of the cluster, can be qualitatively represented in the form Here η = VN p /V cl is a geometric parameter characterizing the packing density of magnetic nanoparticles within the cluster volume, If this hypothesis is correct, the properties of an assembly of non-interacting clusters with a sufficiently large number of particles, N p >> 1, depend mainly on two dimensionless parameters, i.e. the ratio K/M s 2 , which is determined only by the material parameters of nanoparticles, and the filling factor η characterizing the nanoparticle density in the volume of the quasi-spherical cluster.

III. Results and discussion
To check the hypothesis stated above it is necessary to analyze the shape of the averaged hysteresis loops of a dilute assembly of nanoparticle clusters depending on the set of dimensionless parameters (η, K/M s 2 ). To do so, we carried out a detailed calculation of the averaged hysteresis loops of the clusters with different sets of initial parameters K, M s , D, R cl and N p . In the simulations performed the material parameters of the nanoparticles vary in the intervals K = 5×10 4 -10 5 erg/cm 3 , and M s = 300 -800 emu/cm 3 , respectively. Thus, the dimensionless parameter K/M s 2 varies within the range 0.08 ≤ K/M s 2 ≤ 1.1. The diameters of the nanoparticles studied are in the range of D = 20 -40 nm, the number of particles in the cluster is given by N p = 30 -120, the cluster density being η < 0.5.
Random clusters with a given number of particles N p were created in this study as follows. First, we generated dense enough and approximately uniform set of N random points {ρ i } within a spherical volume with the radius R cl . so that |ρ i | ≤ R cl , i = 1,2, ,.. N, N >> N p . Center of the first nanoparticle was placed in the first random point, r 1 = ρ 1 . Then all random points with coordinates |ρ i -r 1 | ≤ 2R were removed from the initial set of the random points. After this operation, any point in the remaining set of random points could serve as a center of the second nanoparticle. For example, one can put simply r 2 = ρ 2 . At the next step one removes all random points whose coordinates satisfy the inequality |ρ i -r 2 | ≤ 2R. This procedure is repeated until all N p nanoparticle centers are placed within the cluster volume. As a result, all random nanoparticle centers lie within a sphere of radius R cl , so that |r i | ≤ R cl , i = 1, 2, .. N p . Furthermore, none of the nanoparticles is in a direct contact with the neighboring nanoparticles.
This algorithm enables one to construct random quasi-spherical clusters of magnetic nanoparticles for moderate values of the cluster filling factor η < 0.5. The orientations of the easy anisotropy axes of uniaxial nanoparticles {e i } were also chosen randomly and independently in the upper hemisphere of the spherical coordinates, 0 < θ < π/2.  Figure 1 shows an example of quasi-spherical random cluster of nanoparticles created using the algorithm described above. Only particles whose centers are located in a layer |z| < R near the equatorial plane z = 0 are shown. The overlap of some particles, visible in figure 1, is only apparent. In fact the particles are at different distances from the plane z = 0, and do not intersect with each other.   Figure 2a shows the hysteresis loop of clusters with different amounts of particles of the same diameter. The volume of the spherical cluster varies inversely to the particle number, V cl ~ 1/N p , to keep the filling factor η unchanged. In figure 2b, conversely, the particle diameters are different, but number of particles in the cluster is maintained constant. The volumes of the clusters are also changed so to keep the η value unchanged. As figures 2a and 2b show, in both cases the shape of the averaged hysteresis loops of the assemblies studied is almost the same. Such insensitivity of the averaged hysteresis loops on the number of particles in the cluster in the limit N p >> 1, or on the particle diameters, have been observed for other filling factors investigated. In figure 3 we compare the averaged hysteresis loops of rare assemblies of clusters with the same filling factor η = 0.4 and the same nanoparticle diameter D = 40 nm, but with different values of the material constants K and M s . For the first assembly with K = 10 5 erg/cm 3 and M s = 500 emu/cm 3 the dimensionless ratio is given by K/M s 2 = 0.4. For the second assembly with K = 5×10 4 erg/cm 3 and M s = 350 emu/cm 3 one has nearly the same value, K/M s 2 = 0.41. The hysteresis loops in figure 3 are normalized to the corresponding anisotropy fields, H a1 = 400 Oe and H a2 = 285.7 Oe, for the first and second assemblies, respectively. As figure 3 shows, due to this normalization the calculated hysteresis loops coincide with sufficient accuracy. Note that under the condition K/M s 2 = const the anisotropy field ratio equals to the ratio of the saturation magnetizations, H a1 /H a2 = M s1 /M s2 . Therefore, the loops in figure 3 will also coincide being normalized on the saturation magnetizations M s1 and M s2 of the nanoparticles, respectively.
Next, we study the dependence of the averaged hysteresis loops on the parameter K/M s 2 at a fixed cluster packing density η. As figure 4a shows, at a fixed packing density η and with increasing K/M s 2 ratio the averaged hysteresis loops of dilute assembly of clusters in variables (M/M s , H/H a ) is gradually approaching the standard Stoner-Wohlfarth [3] hysteresis loop for nonoriented assembly of spherical single-domain nanoparticles with uniaxial anisotropy.   The numerical simulation of the averaged hysteresis loops of dilute assemblies of clusters of magnetic nanoparticles has been performed for the different pairs of dimensionless parameters (η, K/M s 2 ) shown in figure  6. As numerical simulations show, to the left of the curve 1-1 in figure 6 the averaged hysteresis loops are close to the ultimate loops shown in figure 5. The latter are realized in the strong magneto-dipole interaction, M s η/H a >> 1. On the other hand, to the right of the curve 2-2 in figure 6 the hysteresis loops are close to the standard Stoner -Wohlfarth hysteresis loop [3] for non-oriented assembly of non-interacting single-domain nanoparticles with uniaxial magnetic anisotropy.

IV. Conclusions
In this paper, using numerical simulation the averaged hysteresis loops of dilute assemblies of random clusters of uniaxial single-domain magnetic nanoparticles have been investigated. The assemblies studied are characterized by a set of material, K and M s , and geometric, D, R cl , N p parameters. It is shown that the properties of dilute assembly of identical random clusters of nanoparticles are determined by two dimensionless parameters, i.e. the ratio K/M s 2 , and the average particle concentration within the cluster, η = VN p /V cl . In the strong interaction limit, M s η/H a >> 1, the ultimate hysteresis loops of dilute assemblies have been constructed. In the variables (M/Ms, H/Ms) these hysteresis loops depend only on the particle concentration, which varies in the range η < 0.5. On the other hand, in the weak interaction limit, M s η/H a << 1 the assembly hysteresis loops are close to the standard Stoner-Wohlfarth hysteresis loop of non-oriented assembly of uniaxial single-domain particles. Note that some of the restrictions of the model considered in this paper can be safely avoided. Indeed, initially it was assumed that each isolated cluster contains approximately the same number of magnetic nanoparticles N p . But then it was found that the properties of the individual clusters depend actually on the cluster volume fraction η. Therefore, the results obtained are also valid for the case of non-interacting clusters with different numbers of the nanoparticles provided by the clusters have nearly the same volume fraction η. Further, the model assumes that the clusters of the nanoparticles have a quasi-spherical shape. Apparently, the properties of clusters of ellipsoidal shape can be described in a similar way, taking into account the change in the average demagnetizing field existing within the cluster. The model considered in this paper allows generalizations in various directions. It would be desirable to study the changes that occur when one takes into account a distribution of the nanoparticles diameters. It would be interesting also to consider the case when a certain amount of nanoparticles of the cluster can be in direct contact, so that the exchange coupling exists among these nanoparticles. However, these issues deserve special consideration.