Equilibrium properties of assembly of interacting superparamagnetic nanoparticles

The stochastic Landau–Lifshitz equation is used to investigate the relaxation process and equilibrium magnetization of interacting assembly of superparamagnetic nanoparticles (SPMNPs) uniformly distributed in a nonmagnetic matrix. For weakly interacting assembly, the equilibrium magnetization is shown to deviate significantly from the Langevin law at moderate and large magnetic fields under the influence of their magnetic anisotropies. For dense assemblies with noticeable influence of the magneto-dipole interaction, a significant dependence of the initial susceptibility on the assembly density is revealed. The difference between the initial susceptibility and the corresponding Langevin susceptibility can serve as an indication of appreciable influence of the magneto-dipole interaction on the assembly properties. A new self-consistent approach is developed to explain the effect of mutual magneto-dipole interaction on the behavior of dense assembly of SPMNPs. The probability densities of the components of random magnetic field acting on magnetic NPs are calculated at thermodynamic equilibrium. The self-consistent probability densities of these components are found to be close to Gaussian distribution. A decreasing equilibrium assembly magnetization as a function of its density can be explained as a disorienting effect of the random magnetic field on the NPs magnetic moments.

Scientific RepoRtS | (2020) 10:13677 | https://doi.org/10.1038/s41598-020-70711-w www.nature.com/scientificreports/ statistical integral by the Born-Mayer method 55 were employed. However, despite several approaches used, an understanding of the role of magneto-dipole interaction in equilibrium properties of dense SPMNPs assembly is still incomplete. In this paper, the equilibrium magnetization of an assembly of interacting SPMNPs uniformly distributed in a rigid nonmagnetic matrix is calculated by solving the stochastic Landau-Lifshitz (LL) equation [56][57][58][59][60] . This approach is an alternative to the classical method of Gibbs assemblies. It enables one to simultaneously take into account the effect of various types of magnetic anisotropy, magneto-dipole interaction, and thermal fluctuations of the particle magnetic moments on the assembly behavior. Moreover, this method allows one to consider also kinetic processes such as the relaxation process to the equilibrium assembly magnetization.
Calculations based on the stochastic LL equation were performed in this work for a dilute assembly of NPs clusters with a finite filling density η = N p V/V cl , where V = πD 3 /6 is a volume of spherical nanoparticle of diameter D, N p being the number of NPs in a cluster of volume V cl . The random positions of the NPs in the cluster are assumed to be fixed, the rotation of the NPs as a whole is excluded. The easy anisotropy axes of the NPs are randomly oriented. A saturation magnetization of the particles is taken to be M s = 350 emu/cm 3 , which is typical for iron oxide NPs 3,6,10 . The uniaxial magnetic anisotropy constant K 1 is varied from 6 × 10 4 to 1.5 × 10 5 erg/cm 3 . The numerical simulations are carried out at room temperature, T = 300 K. Therefore, the diameter of spherical NPs is restricted to a range D < 25 nm, to ensure 61 that the blocking temperature T b of the largest NPs is well below the room temperature. A significant dependence of the assembly equilibrium magnetization on the intensity of the magneto-dipole interaction inside the clusters has been revealed.
In addition, the statistical properties of random magnetic field acting on magnetic NPs in a dense assembly of SPMNPs have been studied. Following the Lorentz approach 20 it is shown that the random component of magnetic field acting on a typical nanoparticle of the assembly is determined only by the surrounding NPs located inside the Lorentz sphere. The magnetic moment of a typical nanoparticle, and hence the equilibrium magnetization M eq of the assembly, is calculated self-consistently depending on the total magnetic field acting on the particle. The variant of the self-consistent field approximation developed in this work is shown to describe qualitatively correctly the numerical simulation data for the equilibrium assembly magnetization, obtained by means of solution of the stochastic LL equation, over an effective H 0 ≤ 600 Oe range.

Results and discussion
Dilute nanoparticle assembly. As mentioned in the introduction, there are two important contributions that lead to a difference in the reduced equilibrium magnetization of interacting assembly, m eq = M eq (H 0 ,T)/M s , from the Langevin law [19][20][21] where x = M s VH 0 /k B T is the dimensionless Langevin variable and k B is the Boltzmann constant. These are magnetic anisotropy energy and the energy of the magneto-dipole interaction. Let us discuss the influence of these factors on the equilibrium properties of an assembly separately.
Consider first a relatively simple case of a dilute NPs assembly, η → 0, neglecting the influence of magnetodipole interaction. Here, an equilibrium assembly magnetization can be determined evaluating the Gibbs statistical integral [30][31][32] . The corresponding calculations for randomly oriented monodispersed assemblies of SPMNPs are shown in Fig. 1. The magnetic field dependence of equilibrium magnetization of randomly oriented NPs of various diameters is shown in Fig. 1a. Figure 1b shows that the static magnetic susceptibility of the assembly falls down rapidly over H 0 ≤ 200 Oe. As can be seen from the plots, the static susceptibility substantially depends on average D values. Moreover, as Fig. 1c portrays, a second derivative of m eq with respect to H 0 reveals a pronounced minimum in the region of low, H 0 ≤ 50 Oe, magnetic fields, whose position is a function of average D value.
Let us normalize H 0 to the particle anisotropy field, H a = 2K 1 /M s , in a reduced variable h e = H 0 /H a . Then it can be shown 30-32 that a m eq value in a dilute assembly is a universal function of h e that depends only on the reduced height of the particle potential energy barrier, ). However, as Fig. 1d shows in a limit of small H 0 the m eq (H 0 ) curve coincides with the Langevin function, Eq. (1), for all K 1 values. This is a consequence of the fact that in the limit h e → 0 the expansion is valid 30,31 . As a result, a K 1 dependence of m eq in the region of small H 0 disappears. At the same time, according to Fig. 1d for moderate and large H 0 the difference of m eq from the Langevin function is very significant. It follows from Eq. (2) that a static magnetic susceptibility of the assembly in the low-field region does not depend on the K 1 value. Therefore, it is impossible to determine the K 1 value by measuring the static susceptibility of a dilute assembly, dm eq0 /dH 0 , in the limit h e → 0. At the same time, as we will see later, the static magnetic susceptibility of an interacting assembly differs significantly from the Langevin susceptibility. This important fact makes it possible to evaluate the effect of the magneto-dipole interaction on the equilibrium properties of an assembly.
The noticeable influence of the particle magnetic anisotropy energy on the behavior of dilute assembly of monodispersed NPs in the range of moderate and high H 0 fields, and at temperatures not too high with respect to T b was studied in detail both experimentally 23,25 and theoretically [30][31][32] . The area of parameters H 0 and T, where there is a considerable deviation of m eq from the Langevin law was characterized 23,25 as an anisotropic superparamagnetism. Unfortunately, in a number of recent experimental works (see, for example, Refs. [26][27][28][29], Scientific RepoRtS | (2020) 10:13677 | https://doi.org/10.1038/s41598-020-70711-w www.nature.com/scientificreports/ the experimental m eq data are described by a weighted sum of Langevin functions. In this way, the particle size distribution is taken into account, whereas the influence of the magnetic anisotropy energy is completely ignored.
Assembly of dense 3D clusters. As noted in the introduction, a direct application of the Gibbs principle for calculating the equilibrium magnetization of an assembly of interacting NPs is associated with significant mathematical difficulties. To overcome this difficulty, various theoretical methods were used 25, . The most convincing results were obtained by means of Monte-Carlo simulations [33][34][35][36][37][38][39][40][41][42]46,50,53 for assemblies of interacting SPMNPs uniformly distributed in a nonmagnetic media. However, a known drawback of this method is the difficulty in estimating the actual time for evolution of the assembly in a given magnetic field, as individual Monte-Carlo steps do not correspond to real physical time 33 . As an alternative approach to the problem, in the given paper we use direct numerical simulation based on a solution of the stochastic LL equation [56][57][58][59][60] . Numerical calculations of m eq and static susceptibility of a dilute assembly of dense clusters consisting of N p = 60-100 NPs are carried out in a range of applied magnetic fields, H 0 = 0-600 Oe, the cluster filling density being η = 0-0.3. The details of numerical modeling of the kinetic properties of an assembly of magnetic NPs using the stochastic LL equation are described below in the "Methods" section. Figure 2 shows the magnetization relaxation curves of randomly oriented assemblies of magnetic NPs of various D values in a given H 0 field for different initial magnetization states. In the magnetization distribution designated as Z state, at time t = 0 the NPs are magnetized along the applied magnetic field, whereas for the R initial state the magnetic moments of the NPs are randomly oriented in space. Both initial distributions of the particle magnetic moments differ from thermal equilibrium. Figure 2 shows a temporal evolution of the assembly magnetization for t > 0. It is calculated by solving the stochastic LL equation with a sufficiently small numerical time step Δt with respect to characteristic particle precession time T p 60 . To obtain the complete magnetization relaxation curve of an assembly to the equilibrium state a sufficiently large number of the numerical time steps must be taken. The thermodynamic equilibrium is considered to be achieved when the magnetic relaxation curve approaches a constant value, m eq = M eq /M s , and fluctuates around this value with a small dispersion, as shown in the relaxation curves presented in Fig. 2. To obtain statistically reliable results a large number of numerical experiments, N exp = 100-200, is carried out with the same initial conditions. An average magnetization of a dilute  Figure 2a compares the magnetization relaxation curves from a uniformly magnetized state (Z state) in H 0 = 10 Oe for non-interacting and interacting assemblies of NPs of the same diameter D = 21 nm. To obtain the statistically reliable results shown in Fig. 2a, the numerical simulation data were averaged over N exp = 200 independent numerical experiments. In every numerical experiment N = 3 × 10 6 numerical steps were performed with a small time step Δt = 1.26 × 10 -5 μs, assuming a phenomenological damping coefficient κ = 0.5. In Fig. 2a the relaxation curve for a non-interacting assembly, η = 0, can be described by a time dependent exponent with a single relaxation time τ = 0.2 μs. It approaches a steady value, m eq0 = 0.133, which coincides with the reduced equilibrium magnetization of the non-interacting assembly calculated using the Gibbs formula. At the same time, as Fig. 2a shows, the relaxation curve for an assembly of clusters with a filling density η = 0.278 cannot be characterized by a single relaxation time. To approximate this curve at least two exponents with significantly different relaxation times should be used. Nevertheless, as Fig. 2a shows, this curve also approaches a steady value, m eq = 0.056, at a sufficiently long time. It is reasonable to take this value as the equilibrium magnetization of an assembly of clusters with a filling density η = 0.278 in applied magnetic field H 0 = 10 Oe.
As Figure 2a shows, a magneto-dipole interaction leads to a decrease in the magnetization relaxation time at the fast initial stage, followed by a much slower stage of the full establishment of thermodynamic equilibrium, during which the average magnetization of the assembly already changes relatively weakly. Interestingly, the equilibrium magnetization for the assembly of clusters with a noticeable intensity of the magneto-dipole interaction always decreases compared to that of the corresponding assembly of non-interacting NPs.
This conclusion is confirmed by the data in Fig. 2b,c were the magnetization relaxation curves of various assemblies are shown for different initial Z and R states, respectively. As can be seen from Fig. 2b,c, in accordance with the Gibbs principle the equilibrium state of the assembly in a given H 0 field turns out to be the same, regardless of the type of initial magnetization configuration. It is worth mentioning that the Gibbs postulate is not applicable to study a temporal evolution of the assembly magnetization. Fortunately, it can be done numerically by solving the stochastic LL equation. The equilibrium value of the reduced magnetization of the assembly can be obtained by averaging the relaxation curve over a finite interval of times exceeding the characteristic time of magnetic relaxation, t > τ.  Fig. 2b,c. The results so obtained are plotted in Fig. 3. As cluster filling density η increases, resulting in the increase in the magnetodipole interaction intensity inside the clusters, the value of the equilibrium assembly magnetization decreases. In Fig. 3a-c, the magnetic susceptibility of the assembly, dm eq /dH 0 , in the low-field, H 0 → 0, substantially decreases as a function of η values. For a given set of M s and K 1 values for NPs assemblies of D ≤ 21 nm the reduced equilibrium magnetization at room temperature vanishes in the limit H 0 → 0. These assemblies exhibit typical SPM behavior. At the same time, as Fig. 3d shows, for an assembly of NPs of a larger D = 25 nm, there is a remanent magnetization in the limit H 0 → 0. Therefore, the blocking temperature T b of this assembly exceeds the room temperature value. As a result, the true equilibrium state for this assembly is not reached in a finite evolutionary time. Interestingly, in Fig. 3d the remanent magnetization of the assembly decreases with increasing intensity of the magneto-dipole interaction. Figure 4 demonstrates an universal behavior of the equilibrium magnetization curves for assemblies of NPs with a noticeable intensity of the magneto-dipole interaction, η ≥ 0.2. While the equilibrium magnetization curves of assemblies of non-interacting NPs substantially depend on an average D value (see Fig. 1a), these curves in interacting assemblies practically coincide at the same η value. An exception is the magnetization curve in rather large particles, D = 25 nm, with nonzero remanent magnetization.
To explain this effect, one notes that to an order of magnitude the magnetic anisotropy energy of the particle is W a ~ K 1 V, whereas the characteristic energy of the magneto-dipole interaction of the NPs can be estimated as W m ∼ (M s V ) 2 L 3 av , where L av is the average distance between the NPs of the cluster, which can be estimated from the relation L 3 av = V cl N p . Thus, for the characteristic energy of the magneto-dipole interaction one obtains W m ∼ M 2 s V η . Therefore, the energy ratio W a W m ∼ K 1 M 2 s η is independent of the nanoparticle volume being approximately constant for a fixed η value. Figure 5a plots the reduced equilibrium magnetizations of the assemblies of NPs with the same D = 21 nm, but with different magnetic anisotropy constants. It is noteworthy that the equilibrium magnetization of interacting NPs assemblies differs significantly from the Langevin curve. As Fig. 5a shows in a sufficiently dense assemblies with η = 0.278 the influence of particle magnetic anisotropy on the equilibrium magnetization curve is not significant. In particular, the static magnetic susceptibility of the assembly, dm eq /dH 0 , in the limit H 0 → 0 is practically independent of the K 1 value, similar to the case of assemblies of non-interacting NPs (see Fig. 1d). www.nature.com/scientificreports/ However, the static magnetic susceptibility of the interacting assembly is significantly less than the Langevin value, dm eq /dH 0 = M s V/3k B T 30,31 .
To demonstrate clearly the effect of the magneto-dipole interaction on the equilibrium properties of a SPM assembly, it is of interest to study the equilibrium magnetization curves of an assembly of NPs with a negligibly small K 1 ~ 0. As Fig. 5b,c show, the equilibrium magnetization in the case K 1 = 0 approaches the Langevin curve only in the limit η → 0. Note that the magnetic susceptibility of such an assembly in the limit H 0 → 0 substantially depends on its η value. As can be seen from  Self consistent field approximation. The detailed numerical calculations performed above make it possible to quantitatively assess the change in the equilibrium and kinetic properties of the assembly with an increase in the intensity of the magneto-dipole interaction. However, these calculations do not shed light on the physical cause of such changes. It is clear that in the presence of a magneto-dipole interaction, the magnetic field acting on a typical nanoparticle differs from the magnetic field H 0 applied to the assembly, since the magnetic fields of the surrounding NPs also act on this nanoparticle. In dense clusters, at small distances between the NPs, the magnetic fields of the nearest NPs can be very significant. Therefore, a fundamental interest is determining a probability density of such a magnetic field acting on a typical magnetic nanoparticle in the assembly. In recent years, several approaches are proposed [46][47][48][49][50][51][52] to introduce effective magnetic field acting on a typical nanoparticle in a dense SPM assembly. However, it is shown 53 that the expressions suggested for the effective magnetic fields in some cases are hardly consistent with the Monte-Carlo simulation results. In this paper, we develop another approach to evaluate the effect of random magnetic fields acting in a dense nanoparticle assembly.
Let us consider an effectively large spherical assembly, as schematically shown in Fig. 6, which can be characterized by an average magnetization M eq (H 0 ,T) at equilibrium. Let us select around a typical nanoparticle of the assembly a spherical region (Lorentz sphere 20 ) with radius R L much larger than the average distance L av between NPs. Outside Lorentz sphere one can introduce a nearly homogeneous magnetization distribution close to the average assembly magnetization, <M(r)> = M eq . Then, inside the Lorentz sphere, at least near its center, the magnetic field of external magnetic dipoles is almost completely compensated to zero 20 . Therefore, the magnetic field in the center of the Lorentz sphere acting on the reference particle is created by the surrounding NPs located in the Lorentz sphere.
First, let us analyze the probability density of random magnetic field acting on a typical particle of an assembly with a negligibly small magnetic anisotropy constant, K 1 = 0. As Fig. 5b,c show, in such assembly due to the influence of the magneto-dipole interaction a difference arises between the equilibrium magnetization and the Langevin law. Let H = (H x , H y , H z ) be a vector of the random magnetic field in the center of Lorentz sphere created by the NPs located inside it. Without a loss of generality, one can assume that the external magnetic field H 0 is applied along the Z axis of the Cartesian coordinates. Then, the total magnetic field in the center of Lorentz sphere is given by H t = (H x , H y , H z + H 0 ). Let H t be the module of this vector. It is reasonable to assume that in thermodynamic equilibrium the time-average magnetic moment of the reference particle located in the center of Lorentz sphere is where m L (x) is the Langevin function in Eq. (1). It points parallel to vector H t , so Further, let P(H x ,H y ,H z ) be the probability density of a random magnetic field created by surrounding particles in the center of the Lorentz sphere. Then, the average magnetization of the assembly in the direction of the applied field H 0 is given by  Figure 6. A model Lorentz sphere around a reference nanoparticle in an assembly of SPMNPs.
Scientific RepoRtS | (2020) 10:13677 | https://doi.org/10.1038/s41598-020-70711-w www.nature.com/scientificreports/ Thus, to calculate the equilibrium assembly magnetization in the given approximation it is necessary to determine the probability density of random magnetic field in the center of the Lorentz sphere, created by NPs located inside the said sphere.
For given assembly parameters, a self-consistent value P(H x ,H y ,H z ) can be obtained numerically by conducting a sufficient number of numerical experiments with random spherical clusters of a volume V cl , number of particles N p , and a fixed η value. As will be shown below, the partial probability densities P(H x ), P(H y ), and P(H z ) of the random functions H x , H y , and H z are close to the Gaussian distributions. Due to the random nature of the magnetic field H t , which is the sum of a large number of independent contributions of the magnetic fields of individual NPs, there is a relation To find self-consistent probability densities P(H x ), P(H y ) and P(H z ), an appropriate iterative procedure should be performed. At the first stage of this procedure we consider all particles inside the Lorentz sphere to be magnetized strictly parallel to the H 0 field, so that <M x > = 0, <M y > = 0, <M z > = M s . Under this assumption we obtain the empirical probability densities P 1 (H x ), P 1 (H y ) and P 1 (H z ) of the first approximation in the following manner. A sufficiently wide range of magnetic fields, (− H max , H max ), is divided into a large number of equal intervals, ΔH ≪ H max . Then a sufficient number of numerical experiments N exp are performed in random clusters created independently. A random field H = (H x , H y , H z ) in the center of each cluster is calculated and the relative numbers of clusters with components H x , H y , and H z falling into each predefined interval ΔH are determined.
To obtain the partial probability densities of the second approximation, we generate clusters in the volume of the Lorentz sphere, the particle centers being randomly distributed. The magnetization directions of individual NPs are assigned in accordance with the probability density P 1 (H x ,H y ,H z ) = P 1 (H x )P 1 (H y )P 1 (H z ). Namely, the magnetic field H = (H x , H y , H z ) acting on a specific nanoparticle of the cluster is set randomly, in accordance to P 1 (H x ,H y ,H z ) values. Then, the average magnetization of this particle is determined by Eqs. (3) and (4). In this way, we can assign the magnetization of all 'N p − 1' NPs of the cluster and calculate the total magnetic field acting on the test particle. If we repeat this procedure a sufficient number of times, we can determine the empirical probability density in the second approximation, P 2 (H x ,H y ,H z ). These iterations are repeated until successively obtained probability densities, P i (H x ,H y ,H z ), i = 1, 2, … converge to a certain limit. This limiting probability density is used in Eq. (5) to obtain the equilibrium magnetization of the assembly at a given H 0 value.
To obtain the probability density P(H x ,H y ,H z ) with an ~ 1% accuracy, it is enough to carry out only 3-4 iterations of this iterative procedure. The first iteration thus yields partial probability densities P 1 (H x ), P 1 (H y ) and P 1 (H z ), which are very close to the Gaussian distribution, P(H) = exp −H 2 2σ 2 √ 2πσ , with some empirical standard deviations, σ (1) x , σ (1) y and σ (1) z . As a result of the iterative procedure, we obtain series of standard deviations, σ (i) x , σ (i) y and σ (i) z , i = 1, 2, … which quickly converge to some limiting values, σ x , σ y and σ z . Moreover, due to the axial symmetry of the problem an approximate equality σ (i) y is satisfied at each iteration step. As an example, Fig. 7a shows the evolution of the empirical probability densities P 1 (H x )-P 4 (H x ) for the H x component of random magnetic field during four successive stages of the iterative procedure. To obtain empirical probability density, at each stage of the iterative procedure N exp = 10 5 numerical experiments were carried out in which spherical clusters consisting of N p = 60 NPs of diameter D = 21 nm and cluster filling density η = 0.278 were created. To construct the empirical probability densities, the interval of magnetic fields (− 600, 600 Oe) was subdivided into 120 intervals each of 10 Oe. The particle centers inside the cluster volume were randomly distributed using the algorithm described in the "Methods" section. The particle magnetizations were assigned by means of the procedure described above and using Eqs. (3) and (4). As can be seen from Fig. 7a, the successively obtained P (i) (H x ), i = 1-4, values can be described with a reasonable accuracy by the Gaussian distribution. The empirical standard deviations quickly converge to a constant limiting value. The empirical probability densities for the H y and H z components of the random magnetic field H are of the same form. As Fig. 7b shows, for small H 0 values the limiting empirical standard deviations σ x and σ z turn out to be very close each other. As H 0 increases, they  www.nature.com/scientificreports/ bifurcate, but always σ x < σ z . Moreover, σ x = σ y for the transverse components of the random magnetic field due to the axial symmetry of the problem. Figure 8a,b plot so obtained m eq values over H 0 for assemblies with K 1 = 0. Solid lines represent the results of direct numerical calculation using the stochastic LL equation for NPs with diameters D = 17 and 21 nm, respectively, while the dots show the values calculated in the self-consistent approximation developed. The number of NPs in the Lorentz sphere in the latter case was fixed at N p = 60, and only four cycles of the iteration procedure was carried out for every dot. Here, the maximum difference between the results of two calculations does not exceed 15%, which is due to the presence of the correlation effects. Obviously, the dynamics of the magnetic moments of closely located NPs should be strongly correlated, but this fact is not taken into account in the approximation developed. Figure 8c plots m eq values of random assembly of NPs with D = 21 nm calculated for different numbers of NPs in the Lorentz sphere. An increase in the number of NPs in the Lorentz sphere in excess of N p = 60 does not lead to any noticeable change in m eq values.
As Fig. 7b shows, a difference between the self-consistent standard deviations σ x and σ z is usually small in a wide range of H 0 ≤ 500 Oe. Assuming approximately σ x ≈ σ z = σ and performing calculations in a spherical coordinate system with the polar axis parallel to the direction of the applied H 0 field, one can rewrite Eq. (5) as follows where ξ = HH 0 /σ 2 . In the limit H 0 → 0 this integral is estimated to be where σ(0) is the standard deviation at H 0 = 0. For characteristic values of the standard deviation, σ(0) ~ 100 Oe, the Langevin function m L in Eq. (8) changes slowly. Thus, as Eq. (8) shows, with an increase in σ(0) value the initial magnetic susceptibility of the assembly decreases approximately as 1/σ(0). It can be shown that Eq. (8) Obviously, the decrease in the equilibrium assembly magnetization as a function of its density is due to a disorienting effect of the random magnetic field H. Actually, under the influence of random magnetic field the magnetic moments of the NPs on average deviate from the H 0 direction. Similar calculations of the equilibrium assembly magnetization in the self-consistent field approximation were also performed for random assemblies with K 1 > 0 value. Instead of using Eqs. (3) and (4), in this case one has to assign the magnetizations of the NPs within the Lorentz sphere by means of the corresponding Gibbs principle taking into account the K 1 value and the directions of easy anisotropy axes of various NPs in the formulas given in Ref. 32 . Figure 9a,b plot m eq values over H 0 for an assembly of NPs with K 1 = 10 5 erg/cm 3 , M s = 350 emu/cm 3 for the NPs of D = 17 and 21 nm, respectively. The magnetic field dependences of m eq values obtained in the two different methods turn out to be sufficiently close in the entire 0-500 Oe range of the applied H 0 fields.
For completeness, similar calculations of the equilibrium assembly magnetization have been carried out for dilute assemblies of elongated and oblate clusters of magnetic NPs with aspect ratios D z /D = 2.0 and D z /D = 0.5, where D z and D are the longitudinal and transverse diameters of the spheroidal magnetic cluster, respectively. It is shown that for a given H 0 value, the equilibrium assembly magnetization increases for an elongated cluster with an aspect ratio D z /D > 1 and decreases in the opposite case, D z /D < 1, in comparison with the results for a spherical cluster, D z /D = 1.0. These results are explained by the influence of the macroscopic demagnetizing field which acts inside the Lorentz sphere created in elongated or oblate spheroids.

conclusions
An assembly of single-domain magnetic NPs is a complex physical system whose properties are determined by many factors, such as the distribution of NPs in size and shape, the density of the assembly, and the value of the main magnetic parameters of the NPs. Behavior of the assembly depends also on the properties of the medium where the NPs are distributed, beside the applied magnetic field and the temperature. In contrast to classical plasma or quantum gases with Coulomb interaction 14-18 magnetic particles interact via anisotropic magnetodipole forces. Moreover, a SPMNP is characterized by an induced magnetization, which is nearly zero in the absence of magnetic field acting on the particle, contrary to elementary particles whose electric charge is fixed.
In this paper we study the properties of a SPM assembly of monodispersed NPs in a solid nonmagnetic matrix. The calculations performed take into account magnetic anisotropy and the magneto-dipole interaction of particles, but a contact exchange interaction is ignored between the NPs, as they are protected by thin nonmagnetic shells. This model differs significantly from that describing ferrofluids 46,47,50,52,53 . In fluid NPs can rotate as a whole. In addition, they may redistribute to form chains of particles or dense conglomerates 46,47,50,52,53 .
To realize a SPM regime the sample temperature should be higher than the characteristic blocking temperature T b of NPs. It is important that SPMNPs assembly relaxes to a thermodynamically equilibrium over a finite observation time. The fundamental physical quantity of a SPMNPs assembly is the equilibrium magnetization, M eq = M eq (H 0 ,T), which can be easily measured experimentally [23][24][25][26][27][28][29] . Theoretically, this value can be determined on the basis of the Gibbs principle 14-21 as a derivative of the assembly's free energy with respect to an applied field H 0 . However, a direct calculation of the Gibbs statistical integral for an assembly of interacting magnetic NPs involves great mathematical difficulties. In this paper, a new physically adequate method is used for calculating the M eq (H 0 ,T) value of an assembly by solving the stochastic LL equation [56][57][58][59][60] . In contrast to the Monte-Carlo calculations 25,[33][34][35][36][37][38][39][40][41][42]46,50 , the relaxation process to thermodynamic equilibrium in the assembly can be directly observed using the stochastic LL equation. Detailed calculations of the equilibrium magnetization were performed for dilute assemblies of magnetic clusters containing N p = 60-100 NPs of a given diameter. The intensity of the magneto-dipole interaction inside the clusters can be controlled by changing the cluster filling density η.
In an assembly of weakly interacting NPs it is shown that due to the influence of magnetic anisotropy energy, equilibrium magnetization differs significantly from the Langevin law in the range of moderate and large H 0 www.nature.com/scientificreports/ fields. Nevertheless, in sufficiently small H 0 the K 1 dependence of the equilibrium magnetization disappears. In this area the Langevin formula is valid and describes universal behavior of a dilute assembly. For the assemblies of iron oxide NPs studied here the universal behavior is observed over H 0 ≤ 50 Oe. However, for dense assemblies with a noticeable influence of the magneto-dipole interaction a significant dependence of the initial susceptibility on the density is revealed. A difference of the initial susceptibility over the Langevin value serves as a validity of the influence of the magneto-dipole interaction on the assembly properties. In this paper a new approach to describe the influence of random magnetic field acting on NPs in a dense assembly is proposed. In effective field theories [46][47][48][49][50]52 it is assumed that a typical nanoparticle of the assembly is subjected to some self-consistent magnetic field, which takes into account the influence of the magnetic fields of the surrounding NPs. However, in a real assembly each nanoparticle is under the influence of its own local magnetic field which contains a random component. In this paper the probability densities of the components of random magnetic field acting on a typical magnetic nanoparticle are calculated. The self-consistent probability densities of random field components are described by Gaussian distribution. Thus, the standard deviation in the Gaussian distribution becomes an important parameter of the theory. Knowing the probability density of the components of random magnetic field it is possible to calculate the equilibrium magnetization of the assembly in the given approximation as a function of applied H 0 field. It is shown that the approach developed satisfactorily describes the numerical results obtained for the equilibrium M(H 0 ) curve with the help of stochastic LL equation.
The effect of intense magneto-dipole interaction on the properties of an assembly of magnetic NPs is usually explained 25 either by a change in the characteristic height of energy barriers between potential wells of magnetic NPs, or by some collective processes that simultaneously affect the magnetic state of closely spaced magnetic NPs. Based on Eqs. (5) and (6) in this work it is shown that a decrease in the equilibrium magnetization of an interacting assembly as a function of its density can be explained by the disorienting effect of random magnetic field. This leads, on average, to a deviation of the magnetic moments of the NPs from the applied magnetic field direction. In this connection, it is worth noting that the broadening of spectral lines in a high temperature plasma was successfully explained by the action of a random electric microfield, the statistical properties of which are described by Holtsmark 62 or similar 63 distributions.

Stochastic Landau-Lifshitz equation.
Dynamics of a unit magnetization vector α i of a single-domain nanoparticle i of the cluster is determined by the stochastic LL equation [56][57][58][59][60] where γ is the gyromagnetic ratio, κ is phenomenological damping constant, γ 1 = γ/(1 + κ 2 ), H ef ,i is the effective magnetic field and H th,i is the thermal field. The effective magnetic field acting on a separate nanoparticle can be calculated as a derivative of the total cluster energy The total magnetic energy of the cluster W = W a + W Z + W m is a sum of the magneto-crystalline anisotropy energy W a , Zeeman energy W Z of the particles in applied magnetic field H 0 , and the energy W m of mutual magneto-dipole interaction of NPs in the cluster.
For spherical NPs with uniaxial type of magnetic 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 an applied magnetic field H 0 is given by Next, for spherical uniformly magnetized NPs the magnetostatic energy of the cluster can be represented as the energy of the point interacting dipoles located at the particle centers r i within the cluster. Then the magnetodipole interacting energy is where n ij is a unit vector along 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 Scientific RepoRtS | (2020) 10:13677 | https://doi.org/10.1038/s41598-020-70711-w www.nature.com/scientificreports/ The thermal fields, H th,i , i = 1, 2, ... N p , acting on various NPs of the cluster are statistically independent, with the following statistical properties 56 of their components for every nanoparticle here δ αβ is the Kroneker symbol, and δ(t) is the delta function.
The procedure for solving these equations is described in detail in Refs. 57-59 . Random 3D clusters of NPs. In the Monte-Carlo calculations performed to study the SPMNPs assemblies the nanoparticle positions were randomly generated 36,40 on nodes of simple cubic lattices with a certain lattice parameter. This numerical algorithm can hardly be considered as truly random. In particular, it completely prevents the appearance of numerous assembly configurations where certain NPs turn out to be very close to each other, i.e. closer than the lattice parameter chosen. In the present study the 3D clusters consisting of N p identical magnetic NPs with truly random positions were created using numerical algorithm developed in Ref. 12 . First, a dense and approximately uniform set of N random points {ρ i } was created in a sphere of the radius R cl , so that |ρ i | ≤ R cl , i = 1, 2, ... N, for N ≫ N p . The first random point ρ 1 can be chosen as a center of the first nanoparticle of the assembly, r 1 = ρ 1 . Then it is necessary to remove all points with coordinates |ρ i − r 1 | ≤ D from the initial set of the random points. Any random point in the remaining set of points can serve as a center of second nanoparticle of the assembly, for example, r 2 = ρ 2 . Continuing this procedure, one can assign centers to all N p NPs within the cluster volume. Moreover, none of the NPs of the assembly will be in a direct contact with the surrounding particles. This algorithm allows to create random 3D clusters of magnetic NPs with filling densities η < 0.5. The orientations of the easy anisotropy axes {e i }, i = 1, 2, … N p , of NPs in a random 3D clusters are chosen randomly in a sphere.

Data availability
No data-sets were generated or analyzed during the current study.