The competitive effect of non-magnetic defect and films thickness on the ferromagnetic critical temperature in Ising thin-films

In this work, Monte Carlo simulation was employed to investigate the competitive effect of non-magnetic defects and the thickness on the ferromagnetic behavior of Ising spins in a reduced geometry, i.e. thin-films. The magnetic properties were investigated as functions of temperature, defect concentration, and films’ thickness, especially in the ferromagnetic phase transition region. The finite size scaling was performed via the fourth order cumulant of the magnetization to extract the critical temperatures. From the results, the extracted critical temperatures agree well with previous theoretical investigation, where applicable. With increasing concentration of the nonmagnetic defects, the Ising phase-transition-point slightly shifts towards lower temperature, while the increase of films thickness enhances the critical temperature value. Being confirmed by the main-effect-plot analysis, the increase in thickness has much greater influences on the critical temperature than that of the defect concentration, which could be described in term of the average ferromagnetic interaction spin. As the role of the defect is negligence in the range of considered defect concentration (up to ten percent), it therefore suggests that the preparation of ferromagnetic films can be done in normal operating condition where defects usually occur. It may be not economically worth to aim for the perfectly smooth films when the associated application operates at temperatures away from the critical point.


Introduction
In ferromagnetic materials, one type of the defect that usually found under room temperature preparation is the non-magnetic defect, which could be formed from airinduced-vacancy or the contamination of non-magnetic materials. This non-magnetic defect is the main cause to diminish the ferromagnetic magnetization [1,2]. For instance, previous studies found that an important role of the vacancy type defect is to lessen ferromagnetic interaction, which results in the decrease of the critical temperature TC [1]. As the critical temperature specifies the boundary between ferromagnetic and paramagnetic phases (in ferromagnetic materials) [3], ferromagnets could become paramagnets without the change in temperature by introducing high enough non-magnetic defects into the materials. On the other hand, while increasing non-magnetic concentration tends to reduce ferromagnetic characteristic, the increasing of thickness in ferromagnetic thin-films tends to reflect the opposite phenomena. This is as in thicker films, the surface effect becomes reduces, which implies that the average ferromagnetic interaction becomes greater [4,5]. Previous work suggested that when the films is thick enough, the ferromagnetic behavior becomes that of the bulk limit [4]. Therefore, in terms of fundamental understanding, it becomes of great interest how to control the critical behavior at this ferromagnetic to paramagnet phase transition when both the non-magnetic and thickness are taken into account. In addition, in terms of application, this topic is of technological interest as the magnetic system in a reduced geometry is useful for the magnetic recording industries [6,7]. Since, in normal material preparation, the non-magnetic defects commonly occur and randomly distribute in the structure, the properties of the material are then typically different from the ideal defect-free condition. Consequently, any calculation based on the perfect defect-free materials will lead to inappropriate application design. This then emphasizes the importance of non-magnetic defect as well as the reduced geometry to be key factors to control the critical properties in ferromagnetic films.
Nevertheless, the competitive effect of ferromagnetic films' thickness and the non-magnetic defect on the critical temperature has yet to be investigated in details. Therefore, in this study, a more complete picture of the ferromagnetic phase transition in thin-films under the framework of the stochastic Monte Carlo simulation [4,8,9] was aimed as the objective. The study was firstly done by investigating how the overall ferromagnetic properties depend on temperature, films thickness, and non-magnetic defect concentration. Then, the detailed investigation on the critical behavior was carried out to extract the Curie temperature at thermodynamic limit using the finite size scaling analysis [10]. Finally, the results were quantitatively compared with the previous applicable investigation.  8 8 In this work, we considered ferromagnetic thin-films, which usually inherit very strong anisotropy [11,12]. In such the system, when magnetic spins align their directions into anisotropic axis, the over magnetic energy is minimized, so it becomes unlikely for the spins to point into non-anisotropic axis. Therefore, to represent this phenomenon, the infinite anisotropic Ising model was considered with the energy described by the Ising Hamiltonian in the form where <ij> implies that only ferromagnetic interaction between nearest neighboring spins (si and si) is allowed in the sum. In equation (1), the strength of spin-spin interaction is represented by the exchange interaction J. This interaction was also used as unit of energy, i.e. J = 1. As a result, it redefines the unit of the spin as dimensionless, the unit of energy as J, and the unit of temperature as J/kB , where kB is the Boltzmann's constant. Also, as being Ising spin, the spin si (= 1) represents the direction of that spin at site i th . In this work, the spins were arraying on the lattice points of sizes N = L×Ll. The lattice structure was considered with periodic boundary conditions along the in-plane (xy) directions, but free boundary condition along the out-of-plane (z) direction to imitate the films structure. The films thickness l was ranged from 1 to 10 layers. On the inplane direction, somewhat large linear dimension should be chosen to respect the thin-film's aspect ratio.
Consequently, L was ranged from 10l to 10l+60 (in a step of 10), i.e. at least 10 times the size of the thickness l was considered [4].
To introduce the non-magnetic defect into the system, some Ising spins were set inactive to represent the defect. These inactive Ising spins do not occupy the magnetic moment, i.e. si = 0, and give no energy contribution to the Ising Hamiltonian. Further, the number of non-magnetic defects was based on the non-magnetic concentration c, which was varied from 1 to 10 percent of N. As a result, the total number of active Ising spins become N = N(1−c) [1].
Next, to perform magnetic configuration update in Monte Carlo simulation, an initial magnetic configuration was set with all active Ising spins being in the up direction (+1). Then, the Wolff algorithm [13] was used to update the cluster of spins to their opposite direction. This Wolff algorithm has been known as cluster flip algorithm, which helps to minimize the statistical errors arising from correlation time [14]. In the algorithm, a random active spin is chosen as a seed. Then, all neighboring spins to the seed that have the same direction as that of the seed will be linked to form a group of the same spins with a probability of [13] ( ) where T is the simulation temperature. In the group forming, a uniform random number r is drawn from random number generator and the link is successfully set for r  p. After that, the spins that are just added into the group will be used as new seed spins for adding their neighboring using the same Wolff procedure. The group forming stops when no more spin is added into the group.
After that, all spins in the group are assigned their opposite direction (flip), i.e. si → −si.
From the initial state, the system was simulated until reaching equilibration before measuring any observable, i.e. the magnetization per spin and the internal energy Each measurement was taken when the number of flipped spins in the Wolff update exceeded or equaled to N', and this period of time was called Monte Carlo step, i.e. 1 mcs. For each simulation condition, the simulation was run up to 10 4 mcs (or having up to 10 4 magnetization data). Then, the expectation of magnetization per spin where n' was referred to number of measurement used in the time average, and the magnetic susceptibility were calculated. In this work, simulation conditions were thickness l = 1 to 10 (step of 1), linear dimension L = 10l to 10l+60 (step of 10), and temperature T = 1.5 to 4.0 (step of 0.02) J/kB.
For the extraction of the phase transition temperature, the effective Curie temperatures TC(c) for each defect concentration c were located via the fourth-order cumulant UL [10], i.e.  [10], to extract the Curie temperature at thermodynamic limit.

Results and discussion
From the Monte Carlo simulation, the magnetization m and magnetic susceptibility  profiles as a function of temperature were retrieved. For instance, figure 1(a) shows the temperature dependent magnetization results of the ideal defect-free Ising films system for thickness l varying from monolayer to 10 layers. As is seen, the Curie temperature TC, which separates between low and high temperature phases, or between the high magnetization and low magnetization phases, shifts to higher temperatures with increasing the thickness. To describe, according to the Ising Hamiltonian specified in equation (1), the ferromagnetic system arrives it is lowest energy state when all spins point into the same direction, i.e. the ferromagnetic state. However, on increasing the environmental temperature T, due to the laws of thermodynamics, the heat energy flows from the environment to the system until both have the same temperature T. Then, with the energy gained, the spins in the system have more kinetic energy and this causes random thermal fluctuation to the spin arrangement in the system. When the thermal fluctuation is high enough, the ordered state of the ferromagnetic phase is destroyed (high magnetization phase), and the system arrives in the disordered (random) state of the paramagnetic phase (low magnetization phase). Then, when increasing the films thicknesses, the surface effect becomes lessen as there are more spins in the interior of the films with coordination number of 6 (simple cubic case) compared to those on the surface with coordination number of 5. Consequently, with higher average coordination number, the ferromagnetic interaction between a pair of spins (on the average) gets stronger. Therefore, to destroy this stronger ferromagnetic interaction (or the stronger ordered state), it requires higher thermal energy, e.g. from higher environmental temperature, to take to system from the ordered ferromagnetic phase to the random paramagnetic phase. Note that when the ferromagnetic is about to change phase, there occurs large fluctuation of the magnetization value. This level of fluctuation can be characterized by the magnetic susceptibility  shown in figure 1(b) as it represents the variation of the magnetization value, i.e. see equation (5). Therefore, the temperature at the magnetic susceptibility peak can be used to roughly indicate the critical temperature in finite systems, which also agrees with the temperature where the magnetization curve has its maximum gradient. However, these temperatures are not the Curie temperatures, as this word is reserved for the critical temperature of ferromagnetic to paramagnetic phase transition at thermodynamic limit (infinite sized system). To extract the Curie temperature from finite systems, it requires finite size scaling analysis (to be detailed afterwards).
Then, on the inclusion of non-magnetic defects, an interesting phenomenon was evident. Example for bilayer films can be shown in figure 2. With increasing the defect concentration c from c = 0.00 to 0.10, TC was found to shift to lower temperature. This is expected as the more non-magnetic defect brings about the less of the ferromagnetic interaction per spin due to some pair-interaction being destroyed. As previously mentioned, the less of the average ferromagnetic interaction requires less thermal energy for the phase transition from ferromagnetic to paramagnetic phases. Therefore, the temperature where magnetization has maximum gradient (see figure 2(a)) and the temperature where magnetic susceptibility gets its peak (see figure 2(b)) move to lower temperature when increasing the non-magnetic defect concentration.
As the increasing of films thickness l help increasing the critical temperature but the increase of defect concentration c tends to lessen, it is of great interest to investigate how these two parameters concurrently affect the Curie temperature TC of the films system.  However, to extract the Curie temperature for the simulations done in finite sizes, finite size scaling theory is required. This work considered the use of the fourth-order cumulant of the magnetization UL, where at the Curie point, the UL from different L's will have the same value. This cumulant analysis can be shown in figure 3(a). Nevertheless, as can be seen in the figure, if L's are not large enough, there appears some deviation from L U = [10]. It is therefore customary to extrapolate the temperature at the UL's crossing points to infinite limit. This can be done by choosing the smallest L while the other parameters (l and c) are fixed, and set to L. Then, find the crossing temperatures between UL and other UL's (with L > L) and treated as TC(b) where b = L/L. In this study, the crossing temperatures was obtained by interpolation among data. Then, TC(b) was extrapolated to infinite limit (the thermodynamic limit) by performing the linear fit between TC(b) and ln −1 (b). This extraction of the Curie temperature technique can be shown as an example in figure 3(b), which is from the results of 4-layered films. The TC can be marked from the y-intercept value on the TC(b) axis or the coefficient to the zero-degree term in the linear fit, e.g. see the first number in the parenthesis in figure 3(b). As seen, with increasing the non-magnetic defect concentration, the Curie temperature decreases. Consequently, with taking the films thickness into account, the competitive effect of both parameters on the Curie temperature can be shown in figure 4.  As seen from figure 4, the increase in films thickness tend to have the greater effect on Curie point, especially when l is small. This is practical as from l = 1 to l = 2, the number of spins double and the average coordination number changes from 4 to 5, which enhances the average ferromagnetic interaction per spin of 25% increment. However, for thick films, this rate of increment drops, so TC tends to converge to the bulk value at very thick films. On the other hand, the increase of defect concentration tends to linearly remove the active magnetic spins and hence the average ferromagnetic interaction as well as the Curie temperature should also decrease linearly. To confirm this assumption, the main effect plot was perform and shown in figure 5. This plot was generated by keeping only the considered parameter fixed (either l or c) and average all the results from another parameter. As seen in figure 5, it is clear that the defect concentration in the range 1% to 10% has not much of the effect on the Curie temperature, when compared with the change in films thickness. It could be that this range of defect concentration is somewhat below the percolation threshold so the effect on critical behavior is minute [15]. Nevertheless, ten percent or nonmagnetic contamination or vacancy defect is already substantial to be allowed in materials preparation. Therefore, the findings presented in figure (4), can be used as guideline to estimate the tolerance error of the Curie point variation in high anisotropic ferromagnetic materials, when prepared in a normal non-ideal condition.  8 8 Monte Carlo simulations on the Ising thin-films were performed to investigate the competitive effect of non-magnetic defect and films thickness on the ferromagnetic properties, which are the magnetization, the magnetic susceptibility and the Curie temperature. The magnetization strength was found to increase on increasing films thickness and slight decrease with increasing the non-magnetic defect concentration, which is also confirmed by the main effect plot of the Curie point. Therefore, the findings suggest that on preparing ferromagnetic films for magnetic applications, it can be worry-free from not being ideal smooth films usually occurred in normal material preparation/synthesis.