Ferromagnetic order in dipolar systems with anisotropy: application to magnetic nanoparticle supracrystals

Single domain magnetic nanoparticles (MNP) interacting through dipolar interactions (DDI) in addition to the magnetocrystalline energy may present a low temperature ferromagnetic (SFM) or spin glass (SSG) phase according to the underlying structure and the degree of order of the assembly. We study, from Monte Carlo simulations in the framework of the effective one-spin or macrospin models, the case of a monodisperse assembly of single domain MNP fixed on the sites of a perfect lattice with fcc symmetry and randomly distributed easy axes. We limit ourselves to the case of a low anisotropy, namely the onset of the disappearance of the dipolar long-range ferromagnetic (FM) phase obtained in the absence of anisotropy due to the disorder introduced by the latter.


Introduction
The physics of nanoscale magnetic materials is still a very active field of research both in view of the potential applications, especially in nanomedicine [1], and from the fundamental point of view [2,3]. The fundamental aspects are all the more important that magnetic nanoparticles (MNP) can be synthesized in a wide range of size and shapes and their assemblies obtained with different structures: colloidal suspensions, or ferrofluids embedded in non-magnetic materials where one can tune the interparticle interactions through the concentration, or as powders. Moreover, when the size dispersion is sufficiently narrow, MNP can self-organize in long-range ordered 3D supracrystals, namely crystals of nanoparticles [4][5][6], conversely to the discontinuous metal insulator multilayers (DMIM) where the underlying structure is disordered [7][8][9]. The interplay between mutual interactions and the degree of order of the underlying structure is a key factor in understanding the extrinsic properties of single domain MNP assemblies [3]. Indeed, due to frustration effects, according to whether the MNP are arranged or not with a long-range order at high concentration, a ferromagnetic or a spin glass state is expected at a low temperature [3,7] (respectively super-ferromagnetic, SFM and super-spin glass, SSG according to the nanoscale magnetism nomenclature). Hence, one of the purposes of the modelling is to understand under which conditions the SFM or the SSG state can be reached.
The modelling of the magnetic properties of both nanostructured materials and materials including nanoscale particles is a multiscale problem since, on the one hand, the local magnetic structure within MNP at the atomic site scale presents nontrivial features and, on the other hand, the interactions between MNP play an important role. An important simplification in the modelling occurs for particles of the radius under a critical size r sd depending on their chemical composition, typically of a few dozens on nanometers, where they reach a single domain regime avoiding both multidomains and vortex structures (typical values for r sd are 15 nm for Fe, 35 nm for Co, 30 nm for γ-Fe 2 O 3 [2]). In this case, the properties of the MNP ensemble can be described through an effective one-spin model (EOS) where each MNP is characterized by its moment and magnetocrystalline anisotropy energy. However, both the moment and characteristics of the anisotropy term must be considered as effective quantities taking into account some of the atomic scale features. Hence, the effective saturation magnetization M s essentially differs from its bulk value due to crystallinity and/or surface defects, the so-called dead layer, and is expected to be smaller than the latter. As concerns the characteristics of magnetocrystalline anisotropy energy (MAE), the situation is more complex since even its symmetry, either uniaxial or cubic, can differ from that of the bulk material. It is worth mentioning that the non-collinearities of the surface spins can be represented by an additional cubic term in the MAE in the framework of the EOS [10].
The EOS type of approach appears then as a simplifying level of description but, nevertheless, a necessary step to model the interacting MNP assemblies, especially when dipolar interaction are included. In the framework of the EOS models, the theoretical description of MNP assemblies make use of the standard methods of statistical physics since one deals with a set of particles interacting through a pairwise potential in a one-body potential including both the MAE and the Zeeman term due to the applied external magnetic field. For MNP coated with a non-magnetic layer, the interparticle interaction includes an isotropic short-range term and the dipole-dipole interaction (DDI), since the exchange term coming from the direct contact between MNP (referred to as superexchange) vanishes. Hence, one is faced with the models apparented to the dipolar hard or soft spheres, according to the short ranged potential (DHS and DSS, respectively) widely developed in a more general framework of dipolar fluids including molecular polar fluids to which one adds the MAE contribution in the form of a one-body potential. An important distinction is to be made concerning the structure of the system under study according to whether it is a liquid-like suspension, a ferrofluid, or a frozen system which can be an ensemble of MNP in a frozen non-magnetic embedding medium or a powder sample. In this latter case, the thermally activated degrees of freedom are the MNP moments and averages over frozen structure realizations that are done in a second step.
In highly concentrated systems, the nature of the state induced at low temperature by the collective behavior resulting from the DDI strongly depends on the underlying structure and dimensionality due to the anisotropic character of the DDI. For a well ordered system such as a fully occupied perfect lattice free of MAE or with a MAE characterized by aligned easy axes, a long-range ferromagnetic (face centered cubic, fcc, body centered cubic, bcc or body centered tetragonal, bct) or anti-ferromagnetic (simple cubic, sc) ordered state is reached [11][12][13][14][15][16]. Notice that the solid dipolar ferromagnetic ground state presents a bct structure [17], a result in agreement with the solid state phase diagram of the dipolar hard sphere system [18]. Conversely disordered systems present a spin-glass like (SSG) [7,19,20] as has been experimentally evidenced [21][22][23][24][25]. It is worth mentioning that not only the perfect lattices of DHS or DSS [15,26] orient spontaneously at low temperature or high pressure. This is also the case for the dense liquid state [12,[27][28][29][30] where a ferroelectric (or equivalently a FM) nematic state has been reported in the bulk or in slab geometry.
On a perfect lattice, the structural disorder may be introduced either through the MAE contribution with a random distribution of easy axes or through dilution. In the limiting case of an infinite magnitude of the MAE, the DHS on a lattice reduces to the dipolar Ising model, where the exchange coupling constants J i j follow the dipolar energy form for moments aligned on the easy axes. In the perfect aligned dipole (PAD) case, the ordered/disordered transition takes a SG character when reducing the lattice occupation rate [16], while in the random axes dipoles (RAD), the transition is a SG one whatever the lattice occupation [31,32].
Well ordered MNP supracrystals with a fcc structure [4,5] are presumably good candidates for experimental realization of a dipolar induced super-ferromagnetic (SFM) phase. The SFM phase was already looked at through the behavior of the hysteresis in terms of temperature on supracrytals made of iron oxide nanoparticles [6]. A ferrimagnetic phase was evidenced involving the ferromagnetically ordered MNP core moments and the spin canting at the MNP surface. In this work we investigate from Monte Carlo simulations the effect of the magnitude of the MAE on the dipolar SFM state reached by an ensemble of DHS on a fcc lattice when the MAE easy axes are randomly distributed. We also show that the partial alignment of the easy axes leads to the recovery of the SFM state. We only focus on the SPM-SFM transition.

Model and Monte Carlo simulations
We model an assembly of MNP uncoupled exchange, interacting through dipolar interactions (DDI), characterized by their magnetocrystalline anisotropy (MAE) and self-organized on long-range ordered supra crystals of fcc symmetry. Although this simple model is not at all specific, we have in mind more precisely the case of Co or γ-Fe 2 O 3 MNP. The MAE is assumed to be of uniaxial symmetry. This is justified from a crystallographic point of view in the case of hcp-Co, and from the general experimental finding for γ-Fe 2 O 3 where the shape or the surface effect leads to an effective uniaxial anisotropy. The easy axes distribution is random unless otherwise mentioned. In the framework of an effective one-spin model, we consider a monodisperse ensemble of up to N = 2048 dipolar hard spheres located on the sites of a fcc lattice with the volume fraction φ. The total energy in a reduced form, after introducing a reference temperature T 0 and β 0 = 1/(k B T 0 ) reads wheren i andm i are the easy axes and the unit vectors in the direction of the moments, respectively. The coupling constants are related to the saturation magnetization M s , the uniaxial anisotropy constant K u and the MNP volume ǫ d is a usual dipolar constant. We also introduce more relevant coupling parameters, where d nn is the next neighbor distance and a reduced temperature which is a way to introduce the volume fraction, or density, in the system through (d/d nn ) 3 and in well ordered lattices, d = d nn corresponds to the maximum volume fraction (Φ = √ 2π/6 for the fcc lattice). We perform Monte Carlo simulations in the framework of the parallel tempering scheme [33] with a distribution of temperatures bracketing the superparamagnetic (SPM)-ordered phase transition. The set of temperatures {T n } is chosen either from a geometric distribution, T n+1 /T n = const or according to the constant entropy increase scheme [34] which leads to a more homogeneous rate of configuration exchange in the tempering scheme. Periodic boundary conditions and Ewald summations are used with the so-called conducting external conditions, according to which the system is embedded in a medium characterized by an infinite permeability, which eliminates the shape dependent depolarizing effects [13,35]. We emphasize that this is a necessary condition to get a spontaneous magnetization. At this point, we refer a reader to the analysis and comments in [27,36].
We calculate, on the one hand, the mean value of the energy, the spontaneous magnetization, per particle and the nematic order parameter λ [13], revealing the occurrence of an orientational order. The mean value of the magnetization, defined in equation (2.4) is the order parameter of the SFM phase. When an anti-ferromagnetic phase is expected, as is the case for the simple cubic (sc) lattice, the staggered magnetization should be considered in place of m with where p iα denotes the layer in the direction α to which the lattice site i pertains.

33703-3
On the other hand, we calculate the heat capacity C v and the susceptibility from the polarization fluctuation and the location of the C v and δm 2 versus temperature peaks. Here, as in [38], B m corresponds to the normalization such that its limiting values in the (anti)-ferromagnetic and paramagnetic phases are B m = 1.0 and B m = 0.0, respectively. Since we consider a random distribution of easy axes, by increasing λ u we gradually increase the amount of disorder. Starting from λ u = 0 to λ u ≫ λ d , the model goes from the perfect lattice of dipoles to the random axes dipoles (RAD) model where the moments are assigned in the direction to the easy axesm i = s ini and behave as Ising spins s i = ±1. The former presents a well ordered ferromagnetic phase at low temperature [13], while the latter presents a spin glass state (SSG) [31]. According to this picture, we expect the system to present a SPM → SFM and a SPM → SSG transition for small and large values of λ u , respectively.

Results
In what follows, for the sake of simplicity, we shall denote the SPM and SFM as paramagnetic (PM) and ferromagnetic (FM) phases, respectively. Our purpose is twofold: first, the determination of the PM-ordered phase transition temperature in the weak λ u regime, and, secondly, the characterization of the ordered phase. In this work we do not focus on the SSG phase which is known to occur at large values of λ u [31].
In the simulations, we use the parameters of cobalt MNP of ca 8 nm in diameter and the volume fraction corresponding to the experimental supracrystals elaborated in [4] in order to fix the dipolar coupling constant. Doing this, we get λ d ≃ 2.15. On the other hand, the MAE coupling constant is seen as a variable parameter. It is clear that the conclusions hold for systems characterized by other values of the dipolar coupling, since we can adjust the value of the reference temperature T 0 . Taking into consideration the fact that the system free of anisotropy presents a FM phase [13][14][15], we expect the ordered phase to keep its FM character at small values of λ u or at least to present the characteristics of a quasi-long-range order (QLRO) magnetic phase as it is the case for the random anisotropy model (RAM) in the weak anisotropy regime [38].
First of all, in the absence of anisotropy, we determine the PM/FM transition temperature in the framework of the finite size scaling analysis, namely from the crossing point in T * of the Binder cumulant given by equation (2.7), corresponding to different system sizes, see figure 1. We locate the PM/FM transition at T * c = 0.618 ± 0.005. This result, leading to ρµ 2 /(k B T c ) = 2.283 where ρ is the number density, introduced by relating (d/d nn ) 3 of equation (2.3) to ρ, is in close agreement with the early finding of Bouchaud and Zerah [14]. Moreover, as shown in figure 1, we find that T * c is very well bracketed by the C v and χ m peaks, namely T * 02. A precise determination of B m at the PM-ordered transition, namely its non-trivial fixed point value, B * m , is beyond the scope of the present work. It is worth mentioning that this is a difficult task, which is not unambiguous since it depends on non-essential parameters of the transition [26,[39][40][41]. Nevertheless, as expected we get a value, B * m = 0.76 ± 0.02 at λ u = 0, compatible with that of either the Heisenberg model on the lattice (0.7916 [42]) or the dipolar hard sphere model (≃ 0.772 to 0.812 [26]).
Then, we compare the behavior of the simple cubic lattice sc to the behavior of the face centered cubic fcc lattice. It is well known that the sc lattice is antiferromagnetic [11,43], and thus one should compare the staggered spontaneous polarization, given by equation ( of the fcc. The result for both the polarization and the nematic order parameter λ, displayed in figure 2 for simulation boxes of similar size shows that the ordering behaviors of the two lattices present very similar characteristics, at least in the absence of anisotropy. This differs from the situation of the Ising PAD model corresponding to λ u → ∞ with aligned easy axes in the z direction where the PM-AFM and PM-FM transitions are found at T * c ≃ 2.5 and 4.0 for the sc and fcc lattices, respectively [43]. Then, we have performed simulations for the fcc lattice and λ u in the range 0 λ u 1.28 which is still in the low anisotropy regime. We still determine the PM-ordered transition temperature from the finite size behavior of the Binder cumulant and the location of the C v and χ m peaks. An important result is a very weak dependence of the transition temperature with λ u ; indeed, we get T * c = 0.62 ± 0.007 and 0.610 ± 0.007 for λ u = 0.465 and 1.28, respectively. The nature of the ordered phase can be deduced from the low temperature behavior of the magnetization m and the shape of the C v (T * ) curve. We expect the magnetization in the FM phase, on the one hand, to continuously increase with a decrease of T * down to vanishing temperatures and, on the other hand, to be size independent below a threshold temperature. In particular, the continuous increase of m down to T * ≃ 0 is indicative of the absence of freezing of the system at any finite temperature. These two criteria are clearly reached in the absence of anisotropy, as can be seen in figure 3. On the other hand, a lambda-like shape of C v together with the size dependence of the C v peak, as displayed in figure 1 for λ u = 0 is indicative of the onset of QLRO state from the PM one. Indeed, the C v curve shape of typical SG strongly differs from a lambda-like one and, moreover, does not follow the characteristic finite size scaling of the PM/FM (or PM/AFM) transition [16,32]. For λ u = 0.465, we get a behavior for both the magnetization m and the nematic order parameter λ very close to that observed in the absence of anisotropy, λ u = 0 and the above two criteria are satisfied. Moreover, the dependence of m with respect to the number of particles, N is negligible, at least for N 2048, indicating that the mean spontaneous polarization at low temperature remains finite on an infinite system. Therefore, we can conclude that the ordered phase is FM in this case. This is obviously no more the case when λ u = 1.28, as shown in figure 4 where we compare the corresponding spontaneous polarization in terms of T * for different system sizes. Increasing λ u , up to λ u = 1.28, the behavior of m strongly differs and the first point is the appearance of a plateau in the m (T * ) curve at low temperature.  We also get a different behavior for the Binder cumulant, displayed in figure 5. Indeed, if there is still a crossing point at the onset of the PM-ordered phase, (which is, however, less pronounced), the low temperature first B m no more reaches the limiting value B m = 1, and its value decreases with an increase in N which results from the onset of a decrease of B m with T * beyond some value of N. Such a behavior resembles the one obtained in [44] where it was analyzed as the occurrence of a reentrant SG phase. It is worth mentioning that the onset of the PM-ordered transition still presents a FM character, as indicated by the shape of the C v curve displayed in figure 5 but it is presumably a QLRO one instead of a FM phase as indicated first by the size dependence of C v , strongly weakened when compared to the PM-FM transition case. To confirm this point, we show in figure 6   (Color online) Spontaneous polarization (squares) and nematic order parameter (circles) for a system with a textured distribution of easy axes centered on the z-axis with a solid angle limited by a gaussian with variance Θ m = π/6; λ u = 1.86.
the limit N → ∞ for λ u 0.5. From the low temperature moments configurations obtained with 2048 particles and λ u = 1.28, we cannot conclude on the onset of the domain formation at least in this range of sizes; on the other hand, the μ i .μ j correlation function still presents a FM character with no sign inversion but a noticeably reduced range when compared to the λ u = 0 case. Finally, we show that the FM ordered phase is recovered when we consider a preferentially oriented easy axes distribution, with probability p(Θ) ∝ sin(Θ)e −Θ 2 /2σ 2 around the z-axis where one recovers a usual random distribution for σ ≫ 1. In figure 7 we show an example with σ = π/6 in the case λ u = 1.86 a value for which the long-range FM order is lost when the easy axes is randomly distributed.
To conclude, we have shown that an assembly of MNP located on a perfectly ordered fcc lattice presents a FM phase at low temperature only in a reduced range of MAE amplitude, λ u 0.5 for a random distribution of easy axes. However, there is an onset of a QLRO phase, at a temperature T * c very close to the PM-FM transition of the system free of anisotropy. This QLRO phase is characterized by a spontaneous magnetization vanishing in the limit of an infinite system. The FM phase can be recovered by aligning, at least partially, the easy axes.