Grain boundary anisotropy on nano-polycrystalline magnetic thin films

Grain boundaries in polycrystalline thin films with crystallite sizes at nanoscale presents regions characterized by a high degree of local structural disorder. As a consequence, great values of the associated local anisotropies are expected. On this regard, a systematic investigation of the effect of the grain boundary anisotropy on the magnetic properties in such type of nanostructured systems is addressed. For developing this work, a standard Monte Carlo simulation in the framework of classical Heisenberg spins was carried out, with a Hamiltonian involving exchange couplings, dipolar interactions, Zeeman interaction and contributions of cubic magneto-crystalline anisotropy. A quantification of local structural disorder was considered. Results revealed that i) by keeping the same number of grains, different organizations give rise to different spontaneous magnetizations, ii) the critical exponent of the magnetization differs of pure models, which is attributed to the complexity of the lattice and consistent with a distribution of critical temperatures, iii) Boundary anisotropy varies with temperature and its strength are determinant factors for blocking temperatures, and iv) Boundary anisotropy inside in the hysteretic properties where coercive field variations are observed.

New technological components used to confine and guide magnetic fields, such as inductive sensors, flexible antennas and magnetic cores devices, have been developed to take advantage of nano-magnetic properties. The nanostructured systems composed by nanoparticles and ultra-thin films with nano-grains structure had taken great importance by the new and intrigant phenomenally. One of the well-known properties of nanoparticles systems is the fact that coercive field decreases strongly with diminution of the mean diameter in the single-domain regime per particle. Coercivity can goes strictly to zero in a superparamagnetic state (SPM) 1 . However, SPM is not clear in nanograins, the reason is nanograins do not loss completely the connectivity between them. In this case, domain structure can be understood in terms of the magnetic moments fluctuations per nanograin or a low number of them around the superparamagnetic limit (SPL) 2 . Different models have been proposed to quantify the magnetic couplings between grains [3][4][5][6] . The random anisotropy model (RAM) proposed by Herzer is one of the most relevant 5 . This model had contributed to explain some of the properties at micrometric scale. However, new models are necessary at a more reduced scale where the ratio between atoms in grain boundaries and atoms in the core grains increases considerably.
This work presents results of a Monte Carlo simulation of magnetic polycrystalline samples. The study is based on the fact that boundaries alter the magnetic behaviour. The parameters of the different magnetic contributions were adjusted to experimental values. Typically, the magneto-crystalline anisotropy, the surface anisotropy and the dipolar interaction are in µeV range. Therefore, the difference is well stablished respect to exchange interaction which is in the meV range. In case of boundary anisotropy, this might have different orders of magnitude. When both, exchange and boundary anisotropy are in a direct competition, a loss of length correlation in boundaries is reflected in magnetic properties. In particular, the blocking temperature studies give indications of the factors for which the magnetic domains fluctuant independently by thermal effect and the hysteresis loop studies give indications the effect of the intensity of the boundary anisotropy upon the loss of magnetic coupling between grains.

Model description
Polycrystalline thin film samples were simulated by considering a simple cubic (SC) lattice structure. Periodic boundary conditions (PBC) were implemented for giving continuity to the granular growth along the x-y plane, while free boundary conditions along z axis were taken to account. A typical simulated sample is shown in Fig. 1. The procedure for sample construction was detailed in a previous work where a process of structural relaxation in boundaries was taken into account 7 . Linear dimensions are given in magnetic unit cells (muc) according to magnetic moments positions into the crystal lattice. Sample dimension was set at = L m uc 100 in x-y plane with a thickness of = d muc 20 . It was employed 2 × 10 5 atoms in 10 different samples where number of grains ( ) N g was keep at 30.
The code was compiled using gFortran and OpenMP software. A parallelized Monte Carlo method using 10 cores with shared memory access was implemented. In such a way that sample was divided into 10 × 10 × 2 cubic cells of the same volume, each of one was assigned to single computational core. The number of Monte Carlo steps (MCS) was fixed at 3.2 × 10 4 per cell, which was enough for thermalization purposes according to energy relaxation. Magnetization was recorded by following a cooling process where temperature ranged from 400 K down to 2 K every 0.25 K.
The Hamiltonian describing the interactions in the system reads as follows: Where  exc refers to exchange interaction,  an represents the total magneto-crystalline anisotropy which in turn can be broken down in different components,  dip stands for magnetic dipolar interactions, and  h represents the Zeeman interaction due to the influence of a uniforms external magnetic field. Exchange interaction is shown in Eq. 2, where  S i are three-dimensional unitary classical Heisenberg spins and J exc is the corresponding exchange interaction. The sum runs over each i atom interacting with its j th near neighbors within a cutoff radius of muc 3 .
The magnitude of J exc was calculated as a function of the pair distance R ij in the framework of a RKKY approximation [8][9][10][11] . Election is based in our interest of considering variations in ion distance. Such election for J exc is supported by several DFT studies in metals where curve of J exc vs. distance between magnetic ion presents similar tendencies 12,13 . The length of the Fermi wave vector, k F , was set to 1, which is a typical value for metals 11 . J o is a fitting parameter which is chosen depending on the system to be considered. In our case, and for general purposes, this value was fitted in such a way to obtain  www.nature.com/scientificreports www.nature.com/scientificreports/ Components of the crystalline anisotropy term are specified in Eq. 4. Three components were considered: cubic magneto-crystalline anisotropy  cryst , surface anisotropy  surf and inter-granular boundary anisotropy boun  .
It is important to stress that the temperature dependence of the anisotropy has been also considered in attention to experimental works reported [14][15][16] as well as the cubic anisotropy distortion on the surface and grain boundaries, where an important local structural disorder is expected. Studies obtained by ab-initio calculations and experimental processes have shown that the electronic configuration of the atoms belonging to these regions is significantly different respect to the extended and homogeneous crystalline medium in the core of the grains [17][18][19][20] . Different works have concluded that the magnitude of such an interaction is proportional to the number of surface atoms and small local strains may give rise to an important anisotropy contribution 19,21 .
The  cryst term of a cubic nature is presented in Eq. 5, where α i1 , α i2 and α i3 are the director cosines of the magnetic moments respect to the easy axis of magnetization [100] [001] and [010]. Parameters K i 1 and K i 2 account for the magnitude of the mangetocrystalline anisotropy interaction for which a functional dependence on temperature, in addition to a local structural term, has been proposed. Hence, an effective expression is shown in Eq. 6, where the anisotropy values K ni involve the product between an effective temperature dependence K T ( ) n ef 14,22 and an effective crystalline disorder dependence ( ) K v ef i . The sub-index n represents 1 or 2 for K 1 or K 2 respectively and v i is the magnitude of a distortion vector.
On the other hand, anisotropy decreases rapidly as the temperature approaches the critical temperature [14][15][16] . On this basis, we propose the relationship given by Eq. 7 accordingly with experimental reports 11,23 . The parameter T A refers to an inflexion point generally observed in the curves. K no and ∞ K n represent anisotropy values at zero and high temperature. That values were = .
The basis for understanding this surface effect has been exposed by Néel 24 . Where a sum of uniaxial anisotropy contribution per each nearest-neighbor is considered. This approach has been employed for explaining surface effects in nanoparticles 25 and single clusters 19 . Thus, we propose the introduction of ( ) K v ef i presented in Eq. 8 in order to account for the loss or lacking of neighbors (or dangling bonds) and local distortions around each atom i, which become more relevant at grain boundaries and intergranular regions where the degree of local structural disorder is higher and where the crystal symmetry breaks down. Parameter v i is the magnitude of the resulting vector  → R ij of the first neighbors positions which determines the loss of cubic symmetry and it is obtained from Eq. 9.
The parameter γ can be taken as depending on sample properties. In this study we have set γ = 1. The → v i vector plays an important role for surf  and boun  . The magnitude is an indication of the degree of local structural distortion and the resulting direction stands for the single site uniaxial anisotropy direction. A uniaxial surface anisotropy dealing with the surface of the film is suggested in Eq. 10 where an effective surface constant is given by Uniaxial axis is defined by the unitary vector associated to → v i and a parameter of proportionality ε S was fixed in 0.2 according to a previous work 26 , which is a measure of the strength of the anisotropy or the spin-orbit coupling.
In similar way, a uniaxial anisotropy is considered for inter-granular regions. Hamiltonian is represented by The ε B parameter can be fitted according to experimental properties. In fact, the magneto-elastic energy values are considerably larger than the volume magnetocrystalline anisotropy and they exhibit similar magnitudes to the exchange interaction 27 . As a consequence, small strains may give rise to Scientific RepoRtS | (2020) 10:5041 | https://doi.org/10.1038/s41598-020-61979-z www.nature.com/scientificreports www.nature.com/scientificreports/ important magnetic changes 28 . A value of 10 for ε B was used as in previous works 27 , however different values of this parameter can be assumed in order to evaluate the effect of the strength of the anisotropy in the grain boundaries where a high degree of structural distortion is expected.
Concerning the magnetostatic interactions of nano-granular systems, some models have been proposed for explaining the collective behaviour but none of them has been conclusive 29 . Here, dipolar interactions were also considered by using Cartesian coordinates at two levels. The advantage of Cartesian coordinates was demonstrated by comparison to spherical harmonics in Fast Multipole Method (FMM) and Fast Fourier Transform (FFT) in 30 . The size of each cell was the same of the division employed for parallelization purposes. The dipolar energy is given by:   It could be thought that the number of parameters selected for the simulation is high. Nevertheless, all of them are essentials to determine the behaviour of these complex systems. Under a reduced set of parameters, previous studies based upon grains size evidenced a range in which correlation between grains is low and phenomena needs to be better studied 32 . Furthermore, based on previous experiences in the simulation area 7,26,32 , it has been concluded that if parameters are adjusted to experimental ranges, how it was done, an analysis of magnetic behaviour of nano-poly-crystalline films with ferromagnetic magnetic cubic cell can be obtained. That is due to the different magnetic contributions acts in ranges of energy well defined. The differences are stablished respect to exchange interaction. Thereby i) magneto-crystalline anisotropy values are found in the µeV range while exchange interactions values are found in the meV range. Then, with grains in nanoscale, anisotropy depends strongly of the loss of correlation by grain size effect more than magneto-crystalline anisotropy. ii) Experimentally, surface anisotropy generated a perpendicular orientation well stablished at very low thickness, less than 10 muc, then that thickness was stablished like criterium iii) dipolar interaction is proportional to D parameter obtained . All these parameters produce dipolar interactions with D values in µeV range. This is not a parameter with great local variability because depends mainly on local spin. iv) Just grain boundary anisotropy might have different magnitudes orders and alters significantly the magnetic properties by direct competition with exchange interaction when a loss of grain correlation is caused. Both can be found in the same energy scale. Then boundary parameters were selected based in an interval in which boundary exerts notable influence. Figure 3 shows four different results of the magnetization cooling process in different samples by keeping the same number of grains but different realizations (different initial random seed numbers). As can be observed, even though transitions are practically at similar temperatures, the low temperature magnetization can be very different. Such a metastability has been reported to occur in nanostructured ferrite samples of Mg 0.95 Mn 0.05 Fe 2 O 4 , where samples having almost an identical particle size distribution can exhibit different spontaneous magnetization values 33 . In that work, authors concluded that differences are strongly influenced by long range interparticle interactions and by local structural disorder giving rise to different realizations of the grain distribution including different crystallite orientations. In our case, as we will demonstrate, local structural disorder can result in zones of the sample, e.g. those with an average coordination number greater than the nominal one, making the short range magnetic coupling between grains a determinant factor for having different spontaneous magnetizations.

Results
Zero-field-cooled (ZFC) and field-cooled (FC) curves for different external field values are shown in Fig. 4a. Each curve is the average over 5 simulations. The sample is first cooled without-field until the lowest temperature. After that, an external field is applied while temperature increase, and the magnetization is recorded. This record is denoted as ZFC. Subsequently, the sample reached the highest temperature, data are recorded during cool process to lowest temperature keeping the external field. This record is denoted as FC. As expected, ZFC-FC curves are smoother than those presented in Fig. 3, due to the preferential anisotropy direction of each grain imposed by the external field direction.
The critical exponent β associated to the magnetization and the transition temperature T C were inferred by fitting the data in a vicinity around T C inspired by concept of critical temperature distribution per grain presented by Berger, A 34 . and the following relationship 35 : where = − t T T 1 / C and A, B, C, β ′ and β″ are fitting parameters. Results are summarized in Figs. 4b,c respectively as a function of external field. Extrapolation to zero field allowed to obtain β = . ± . 0 46 0 03 and = . ± . T 327 3 2 06 C . Our exponent is greater than the one observed in pure and homogeneous 3D classical Heisenberg models having β = . 0 36 35 . Such a difference is attributed to local structural characteristics of our system not observed in pure models (e.g. single-crystal films) and consistent with a distribution of critical temperatures. Figure 4d presents the results of blocking and irreversibility temperatures, T B and T I respectively. The difference between these temperatures is ascribed to the grain size distribution having different T B 36 . Both curves decrease in a monotonous manner. The behaviour found is in agreement with different experimental results presented by M. Knobel et al. 37 and other reports 38,39 . In the former, authors compare the field dependence of the blocking temperature in a γ-Fe 2 O 3 monolayer sample and diluted nanoparticles. While for non-interacting nanoparticles a linear behaviour with  h was observed, for interacting particles the experimental behaviour showed a similar trend to the obtained in this work.
On the other hand, changes in crystalline anisotropy were implemented to analyse the influence over the blocking temperature. Different magnitudes of K 1 did not show any influence. However, changes of T A , the temperature at the inflexion point of K T ( ) n ef (see Fig. 5a) evidenced an interesting behavior as can be observed in Fig. 5b through the ZFC-FC curves. Initial configurations were obtained from previous cooling processes. As it was already pointed out, in the cooling process spontaneous magnetization can reach different values. Therefore, ZFC curves may eventually overlap at some points. Nevertheless, as it is observed in Fig. 5c, the blocking and irreversibility temperatures exhibit a defined tendency to decrease as a function of T 1/ A . Figure 6 shows the correlation between the structural landscape and the local S z components by means of colormaps representations of a final state in a cooling process without external field and at different values of ε B . Such values can vary in order to evaluate the effect of the strength of the anisotropy in the boundaries. To do so, several values of the ε B parameter were considered. At low ε B values, the magnetization is more homogenous without sharp changes of the local S z components as is presented in Fig. 6b. As ε B increases, a clustering magnetic effect appears as it is possible to observe in Fig. 6c for ε B = 20, where such a domain may involve more than one grain. In contrast and according to Fig. 6d, high values of ε B lead to a boundary disorder effect with an inhomogeneous magnetization. Only big grains can present some reduced domain in their interior. A similar behaviour was observed when considering the other components S x and S y but with a greater number of domains involved. This is due to the preferential x-y given by dipolar interaction and the symmetry of the system. Figures 7a,b show the effect of the boundary anisotropy parameter ε B upon the shape of the ZFC-FC curves respectively. These curves are the average over five simulations. The respective blocking and irreversibility temperatures are presented in Fig. 7c as a fuction of the ε B parameter. An increment in both quantities is observed. When increments in the boundary anisotropy are introduced, domains turn out magnetically harder giving rise to ha higher barrier to be overcome and therefore a greater blocking temperature is needed. It is interesting to remark however, that contrary to the well-known linear relationship of the blocking temperature with the density of anisotropic energy for a given volume, here, two different linear regimes seem to be observed, one below ε = 20 B and the other one at higher values. According to the results presented in Fig. 6, such value of ε B corresponds to the limit above which grain boundaries coincide with the domain limits. The m oFC , the magnetization at 2 K in field cooling is presented as a function of ε B parameter in Fig. 7d. According to this figure, high values of ε B make that the disorder effect produced at the borders destroys homogeneous magnetization within the grains. Therefore T B , T I and m oFC response to this variations of ε B , where at low values of ε B the pinning of the magnetic moments in the domains depends on the collective behaviour of the film, at intermediate values of ε B the insulating effect of the borders makes the domains become independent in each grain, and at a high value of ε B that the film behaves similar to a pin glass system, in which the coupling between the magnetic moments shows a strong degree of frustration.
Regarding the hysteretic properties, hysteresis loops were also simulated and they are shown in Fig. 8a for different values of ε B at a temperature below the critical temperature. As can be observed M-H loops the system becomes magnetically harder as the parameter ε B increases, which in turn makes the coercive force to increase as it is shown in Fig. 8b, whereas the remanence tends to slightly diminish (see Fig. 8c). Nevertheless, it is worth www.nature.com/scientificreports www.nature.com/scientificreports/ noting that despite of varying the boundary anisotropy, no humps in the hysteresis loops are observed, which means that the mechanism for magnetization reversal takes place in a gradual way where the system behaves as a whole and not in a differentiated fashion implying separated contributions of the grain cores and grain boundaries. Different hysteresis loop simulations varying T A and K 1 do not exert influence upon the coercive force and remanence. That is due to exchange interaction which is 200 times magnitude orders greater than crystalline anisotropy.

conclusions
The interplay of the temperature dependence of the grain boundary anisotropy and local structural disorder in nanostructured thin films was analyzed. Results revealed that i) by keeping the same number of grains, different realizations gave rise to different spontaneous magnetizations, ii) the critical exponent of the magnetization was different from that of pure models. The difference was attributed to the complexity of the lattice structure in accordance with the distribution of critical temperatures found in other reports of inhomogeneous films iii) the way in which the boundary anisotropy varies with temperature and its strength were determinant factors for blocking temperatures, and iv) hysteresis loops below critical temperature were characterized by a high degree of symmetry with a coherent mechanism of reversal rotation, without humps or jumps.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. OriginPro 8 software was used for image processing 40 .