Spin-glass-like behaviour of a collection of nano-granular magnetic particles embedded in an insulating thin film and interacting only via dipolar forces

We perform Metropolis Monte Carlo simulations of the behaviour of a film of insulating material containing a distribution of magnetic nanoparticles. We assume that these particles only interact through dipolar forces and we find that their behaviour at low temperatures shows characteristics of a spin-glass with a freezing temperature at which the linear susceptibility and the specific heat exhibit a maximum. We define an order parameter which shows a similar dependence on the temperature to the one expected for a spin-glass . We also calculate the correlation of the orientation of the magnetic dipole at the center of the system at different times. All our results are consistent with the temperature dependence of the variance and the mean of the local field acting on the central dipole. Furthermore, the total magnetic dipole of the system scales with the square root of the number of particles, which indicates that the dipole orientations are randomly distributed, as expected for a spin-glass.


Introduction
At present there is a growing interest in low dimensional magnetic systems. In particular, many studies are being conducted on granular systems which consist of nanometric clusters of a magnetic metal dispersed in a nonmagnetic solid matrix. Such systems have received much attention because of their potential application to ultra-high magnetic storage capacity [1]. When the matrix is electrically conducting (metal or semiconductor) the composite system can exhibit collective behaviour due to indirect exchange or Rudermann−Kittel−Kasuya −Yosida (RKKY) interactions [2]. In some cases these systems are ferromagnetic with a critical T c in the room temperature range [3]. Some of these films also display giant magnetoresistance [4].
When instead the host matrix is an insulator, or in general when the charge carrier concentration in the host is very low, the RKKY interactions between the magnetic particles can be neglected. If the average inter-particle separation is larger than a few Ångstroms we can also neglect the effect of both direct exchange and super exchange, since both require that the magnetic particles be in atomic contact [5]. This situation is not included in our simulations. In such a case, which is exemplified by some of the composite films mentioned above, the collective behaviour of the system is controlled by the magnetostatic (dipolar) interactions at low enough temperatures.
Detailed numerical calculations of the RKKY and the dipolar interactions between clusters of varying sizes show that for realistic cases the latter can be as important or even more important than the RKKY indirect exchange interactions for realistic situations and, furthermore, that the magnetostatic interactions can be well approximated by substituting each cluster by a point magnetic dipole located at its center [6].
Experimental work on the dynamics of several different systems in which Fe particles are dispersed shows that dipolar interactions can control the behaviour of the dynamic susceptibilityfor an adequate range of particle Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. diameter and concentration. When the strength of the inter-particle interactions is increased, either by increasing the concentration of the particles or their radius, and hence their magnetic moment, the dynamic behaviour transitions from an interaction-modified super-paramagnetism to a glassy-type collective dynamics [7][8][9].
A simulation of a planar triangular lattice of nanoparticles with random anisotropy and interacting through dipolar forces has shown the interplay of both energies [10].
In this article we use numerical simulations to study an assembly of magnetic point dipoles which only interact through dipolar forces. Each dipole is located at a fixed but random position inside a non-magnetic insulator film, with the view to finding out whether this system exhibits spin-glass properties at low temperatures. The only restriction on the space configurations is that the distance between the centres of each pair of particles be greater than their diameter. We assume that their magnetic moments can orientate freely because we neglect the magnetic anisotropy, since the emphasis is on the effects of the long range dipolar interactions. Each dipole represents a ferromagnetic nano-particle. In order to obtain meaningful values of the physical properties calculated in the simulation one must stipulate the magnetic moment of the ferromagnetic particles which are assumed to be identical. As a realistic example we assume that each ferromagnetic particle exhibits the total magnetic moment of a body centered cubic Fe [11] spherical nanoparticle, which also requires us to define the value of its diameter. We present the results of Metropolis Monte Carlo (MMC) [12] simulations for this system. The model and algorithm are described in the next section and the simulation data are exhibited and discussed in the subsequent sections. Emphasis is placed on the spin glass nature of the results and to this purpose we examined both global properties, such as the the energy, heat capacity, magnetic susceptibility and spin-glass-like order parameter, and local properties such as the auto-correlation function of the central dipole and the local magnetic field induced by the interactions between dipoles.
It is important to note that the simulation results reported in this article refer only to the dipolar contributions to the various physical quantities.We also note that the ferromagnetic nanoparticles could also be composed of other magnetic atoms such as Co or Ni

Model and Algorithm
We incorporate periodic boundary conditions in the film by introducing a square lattice of cells and we then randomly distribute = N n n x z 2 particles inside each cell. Here n y =n x . This basic cell is then repeated indefinitely along the (x,y) plane.
We use in the following bold face typesetting for vectors and tensors when there is no possibility of confusion.
If we take m (i) to be the dipole moment of nanoparticle i, the total dipolar energy, W, of the basic cell of the film is given by: where r n is a site in the square (x, y) lattice and the components of the dipolar tensor, D, are The previous equation defines the local field, B (i) , acting on nanoparticle i, which will also be calculated. The dipole moments of the magnetic particles can be expressed as: where μ 0 is the magnetic moment of each nanoparticle, and u (i) is the unit vector in the direction of the dipole moment of the ith nanoparticle. We denote each particle by a superscript and the components of vectors or tensors by subscripts.
In this model, frustration, which is a necessary ingredient for glassy behaviour [13], is a consequence of the disorder in the positions and orientations of the dipoles, resulting in random sign and amplitude fluctuations of the tensorial dipolar interaction and the local field. We have chosen a fixed value of n z =3 , which corresponds to a width of the film equal to 3×f×d ,where d=10 Å is taken to be the particle diameter and f>1 is a factor which determines the average inter-particle distance. We choose f=2. This is consistent with reported experimental values of average inter-particle distance in similar films [7].
Clearly, for a given spatial distribution of the magnetic particles the energy, W, only depends on the orientations of the N unit vectors u (i) .
The dipolar sum in equation (1) is calculated by the adaptation of Ewald's summation algorithm to the quasitwo-dimensional case [14]. We use the MMC algorithm [12], as described in the appendix, to calculate both thermodynamic and local properties of the system and we finally average all properties over 100-200 random spatial distributions of the particles.
In the following equations, external angular brackets indicate configuration averages and internal ones indicate thermal averages. We place a particle at the geometric center of the repeating cell in all configurations and we calculate the average and the variance of the components of the local field, B (c) , and of the dipole, m (c) , at this site. We assign the number 1 to this site in all simulations. The averages of all properties are extrapolated to  ¥ n n , x y to obtain the limit of an infinite film of width n z =3. Thermal and configuration averages of the following physical variables are calculated as functions of temperature, T, for different numbers of particles, N, in the basic cell: (1) The averaged total dipolar energy of the cell (the internal energy) is given by: where W is defined in equation (1).
(2) The heat capacity is given by: The total dipole moment of each cell is given by: The components of the susceptibility tensor are defined as follows: The averaged time auto-correlation function, g(t), of the central dipole orientation vector is defined as: where the subscript t 0 , indicates that an average over this initial time is performed in addition to the thermal and the configuration averages. Note that in our calculations time t is defined as the number of MMC steps (see appendix) as the simulation progresses. This implies from the definition above that we calculate the scalar product of u (c) , after t 0 MMC steps with its value after a further t MMC steps, as defined in the appendix. The time scale is not otherwise defined, being dependent on the average physical time a dipole takes to change its orientation, which depends on each real system. Furthermore, the time t above is measured in units of the number of MMC steps, so that the decay time of the correlation is also measured in these units. We emphasize that quantitative information is not relevant as our only aim is to characterize the behaviour of the system at low T. From the original Anderson definition of the order parameter q EA [15] we have: EA t An alternative order parameter has been defined as follows [16]: Here we take t w =N w , the number of preparation steps in a warming cycle, and τ=N MC , the number of steps in the calculation proper. The indicated average is over spatial configurations, while the succesive values of u ( i) α are those obtained in the MMC algorithm. We calculate both order parameters for every value of the temperature T.
In what follows we express temperature T in degrees Kelvin (K).

Thermodynamic properties
At T=0 both order parameters defined in the previous section should equal unity while they are expected to vanish above a freezing temperature T f . However due to the system's finite size, both q EA (T) and q H (T) have a long tail for high T as can be seen in figures 1 and 2.
One can fit the curves for q EA and q H versus T with an algebraic function such that they intersect the temperature axis at an extrapolated temperature a ( ) The freezing temperatures obtained from both order parameters are in good agreement with each other, indicating that both definitions are consistent. However, q EA is amenable to a better extrapolation, since it is less sensitive to the finiteness of the sample used in the simulation. At any rate, one has to cut the data off at a  reasonable value of T≈1.4-1.6 K, since, as mentioned above, the order parameters for a finite system do not vanish at any finite T. This arbitrary procedure introduces an uncertainty in T q of the order of 0.1-0.2 K. In figure 3  It is noteworthy that the mean field solution of the EA model of a spin glass, as mentioned by Mydosh [17], yields at low T : T T 1.0 0.406 6 13 EA f in close agreement with our result in equation (12).
We also obtain an estimate of T f by calculating, following Binder [18], the kurtosis of the distribution of the total dipole moment, M, as a function of T, defined as  The temperature at which the plots of this quantity versus T intersect for different values of N should give an estimate of the freezing temperature for infinite size [18]. It should be noted that the kurtosis is zero for a Gaussian distribution.
The average of the total energy of the cell (equation (3) ) is shown in figure 4 for three different values of N. We verify that it converges as N increases. The heat capacity is shown in figure 5 as defined in equation (4), in dimensionless units. Just as for the parallel static susceptibility, χ xx (T), the curve in figure 6 has a maximum at a temperature, T m (N), slightly lower than T q (N). These results are compatible with spin-glass behaviour [19].
The graphs of the Binder cumulent (kurtosis), defined in equation (14), as a function of T for three different system sizes intersect approximately at a temperature, T c , which is considered to be an estimate of T f . From figure 7 we obtain T c ≈1.5 K. We verify (figure 8) that the average modulus of the total dipole moment of the basic cell, defined in equation (5) scales as N 1/2 which corresponds to a random dipole distribution. In order to see why this is so we calculate the square of the sum of unit vectors in equation (5) which can be written as  If we assume that all sites are equivalent, which is true in a very large system, we can write the second term on the right hand side of equation (  The last equation is obtained by comparing this expression with equation (5) from which it can be seen that figure 8 is the plot of + D ( ( )) ( ) T 1 1 2 : at high T the correlations contained in Δ(T) tend to vanish due to thermal fluctuations, while they grow when T is less than T f .

Local properties
The time auto-correlation function, g(t), defined in equation (7) is shown for N=108, T=2.9 K and 160 configurations in figure 9. By definition g(0)=1. The best fit to g(t) was found to be where, as indicated above, both parameters are functions of T.  (17) is shown as a function of temperature in figure 11.  We see that as  T 0 g(t) tends to one while the decay time of the transient vanishes. As T increases above T f , the transient gets longer and the limiting amplitude decreases, and should eventually vanish. We again find an abrupt change in the derivative of the curve for τ(T) which is the signal of the freezing temperature.
As a necessary step in the calculation of the energy and the specific heat we need to obtain the local field, B (c) , on each particle as defined in equation (1). The statistical distributions of the values of all three components of this field are obtained as the MMC run proceeds, and we use them to calculate its average and variance. Figure 12 shows that the modulus of the local field on the central dipole increases as T is lowered, while its variance, shown in figure 13, decreases rapidly for T below the freezing temperature T f ≈1.5-1.6 K The central dipole exhibits a similar behaviour: figure 14 shows that the variance of (u (c) ) x decreases abruptly below T f .
We conclude that below T f each dipole orientates in a given average direction, around which it fluctuates with a decreasing variance as T decreases. This is shown by the change in the derivative shown in figure 14.
On the other hand, both the heat capacity ( Figure 5) and the static susceptibility ( figure 6) show a maximum at about the same T m .

Conclusions
We present a simulation of the behaviour of a collection of nano-particles sustaining a magnetic dipole moment dispersed randomly in a non-magnetic insulating film (the same conclusions are valid for electric dipoles for  instance in a liquid crystal). We focused our study on the range of high concentrations, where collective behaviour can be expected.
We find that such a system exhibits a freezing transition at low temperatures similar to that of a spin-glass, as shown by the temperature dependence of both local and global statistical properties, namely: (a) One can define an order parameter which is unity at very low T and decreases as T increases; (b) Each magnetic dipole at low T tends to orientate in a fixed direction, around which the amplitude of its fluctuations decreases as T is lowered. This is shown by the correspondingly decreasing variance and increasing average amplitude of a given component. Furthermore, the directions of the dipoles are random at temperatures above the freezing temperature, T f , while they become increasingly correlated below T f , which is shown by the fact that the total magnetic moment of the basic cell per particle scales as N −1/2 above T f ; (c) Both the magnetic specific heat and the longitudinal static susceptibility show maxima at about the same temperature, which is close to the estimated freezing temperature obtained by the Binder criterium based on the temperature dependence of the kurtosis of the statistical distribution of the total dipole moment; (d) The local magnetic field on a given dipole exhibits increasing amplitude and decreasing variance as T is lowered below a given temperature which coincides with the temperature at which the same changes occur in the plots of the central dipole components, and which we interpret as the temperature, T f , of the freezing transition.