Micro mean-field Monte Carlo (MMFMC) method for Heisenberg ferromagnet with arbitrary spin and single-ion anisotropy

We proposed a micro mean-field Monte Carlo (MMFMC) method based on Metropolis single spin flip algorithm to calculate the thermodynamic properties of three-dimensional ferromagnet describerd by Heisenberg model. In MMFMC method, the actions of the neighboring spins on the selected spin is replaced approximately by their expected values multiplied by the exchange interactions between them, thus each spin is in a quantum state obtained by solving single spin Hamitonian in every MC step. The magnetization, internal energy, specific heat, longitudinal and transverse susceptibilities are calculated for SC, BCC and FCC ferromagnets with arbitray spin quantum number and single-ion anisotropy. it is found that the system shows phase transition phenomenon from paramagnetic to ferromagnetic phase with temperature decreasing, and the Curie temperatures obtained are consistent with those obtained by other theories.


Introduction
The ground state of three-dimensional (3D)Heisenberg ferromagnet is ferromagnetic state with spin paralleling each other. With temperature increasing to a definite point, the system transitions to paramagnetic state. In paramagnetic phase, the magnetic susceptibility of system follows Curie-Weiss law which was interpreted by Weiss molecular theory. It was not until the discovery of exchange interaction by Heisenberg that the nature of molecular field was clarified, and the molecular field theory was found equivalent to the mean-field theory (MFT) of Heisenberg spin system [1]. After that, the thermodynamics and critical behavior have become one of the key problems of ferromagnetic system, and many research methods are proposed, such as clusterapproximation [2], Green function method(GFM) [3], high-temperature-expansion technique(THE) [4], Monte Carlo (MC) simulation [5] and so on.
The Weiss molecular theory was the first theory to interpret the ferromagnetism of ferrommagnet, in which the interaction of a spin with the other spins around is replaced by a mean-field, where h in is the internal mean-field independent of the site i, the transition temperature, = + T zJS S 1 3, c ( )/ and the critical behavior calculated can be compared with the experimental results. Considering the correlation of spin itself strictly, Kaneyoshi proposed an effective-field theory with correlations (EFT) by using operator equality technique to treat the Ising model in transverse field [6], improved the results of the MFT. Johnston extened the MFT to the noncolinear Heisenberg antiferromagnet [7] to interpret the susceptibility data. It is found that the increase of the number of spins in a cluster can improve the accuracy of calculating physical quantities [2,8]. The Green function method(GFM) proposed by Bogolyubov and Tyablikov for spin-1/2 ferromagnet [9] and extended to arbitrary spin S is suitable for use in the whole temperature range [10,11]. Improving the treatment of the single-ion anisotropic term and decoupling approximation to the exchange term, this method has become one of the effective means to study Heisenberg system for arbitrary spin and single-ion anisotropic strength [12][13][14][15][16][17].
At high and moderate temperatures the high-temperature expansion (HTE) technique is a very general approach to calculate thermodynamic quantities of Heisenberg system [4,[18][19][20][21]. Although the calculation process is very complicated, the critical point and critical behavior obtained are almost the most accurate. With Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the development of computing technology and equipment, computer simulation has become an essential tool to study of phase transitions and critical phenomena for Heisenberg system. Based either on the Suzuki-Trotter approach [22] or on the Handscomb method [5], quantum Monte Carlo (QMC) gave the critical point of Heisenberg ferromagnet very closure to those of THE [19][20][21][23][24][25] for S=1/2. Although it is a promising method, QMC simulations at finite temperatures are considerably more time consuming than the corresponding simulations of classical systems. In addition, this theory is profound limited to the study of systems with large anisotropy and quantum spin system [26].
In traditional MC simulation, the Metropolis sampling algorithm in single spin flip is successful in simulate the classical model, and the concept of the MFT to treat Heisenberg system is very clear. Combined with the above two methods we propose a new numerical method, micro mean-field Monte Carlo (MMFMC) method, to simulate the Heisenberg ferromagnet. In every MC step, the interactions of a spin with other spin around it is replaced by their expected values multiplied by the exchange interaction between them, thus the quantum state of this spin can be obtained easily. In this way, the quantum states of the spins of systm in every Monte Carlo step can be obtained and finally the quantities of the system can be obtained. In this article we calculate thermodynamic properties and Curie temperature of the 3D Heisenberg ferromagnet with SC, BCC and FCC structure. In section 2, we explain the MMFMC with spin model. The calculation results for thermodynamic properties and Curie temperature for several spin quantum number are arranged in section 3. Summary is given in section 4.

Model and simulation method
The Hamiltonian of a ferromagnet is given by where J ij represents the exchange interaction between spins S i and S , j 〈ij〉 represents the summation over the nearest-neighbor spin pair, D i is the single-ion anisotropic constant, h is the external magnetic field applied along the z direction.
The Hamiltonian of spin S i in the system can be written as, In traditional MFT, the effect of the surrounding spin on spin S i is replaced by their thermodynamic mean value mulitiplied by the exchange interaction between them, that is, and the quantity á ñ S can be solved self consistently. So the MFT omits the thermal fluctuation of spins, and the critical temperature obtained is higher than those obtained by the other method. Now if we replace the thermodynamic mean value á ñ S j with the mean value of quantum mechanics á ñ S j of S , j that is, here á ñ S j is dependent on the site j and is a constant if the quantum state of spin S j is given. Thus the Hamiltonian and we can solve this Hamiltonian and get the quantum state and expected value of spin S i , which are the functions of the expected values of the nearest neighbor spins. By using MC simulation based on Metropolis single spin flip, we can get the thermodynamic quantities of the system, we called this method as the micro mean-field Monte Carlo (MMFMC) method. In executing MMFMC method, we first assume that the spins of the system are in their initial quantum states and system is in the superposition state of these states, and the quantum expected value of each spin in the lattice, á ñ S , i i=1, 2,K, N, can be determined in this quantum state, and so are the energies of each spin = á ñ = E H i N , 1, 2,..., .
i i Then we get the eigen-state and the eigen-enegy of the spin by solving the single spin Hamiltonian (4), and accordingly we get the expected values á ñ S i n of the spin S i in these states with á ñ = á ñ = + n n n S S S , 1,...,2 1.
i n i | | By comparing the directions of the vectors á ñ S i n and á ñ S , i we exclude the state closest to the initial state, and select randomly a state in the remaining 2 S states to compare with the initial state, we get the new quantum state of S i based on Metropolis algorithm. If á ñ » S 0 i in the initial state for system with spin quantum number S being even number, the state closest to the initial state in the 2 S+1 states is the state in which the expected value of S i is approximately 0. Finally the new state of system is obtained. To go back and forth, we get a series of states of system considering the thermodynamic weight. For Ising model, because the commutation between spins, the spin Hamiltonian is diagonal, the MMFMC method is changed to the traditional MC method.
In MMFMC method, although the eigen-states and expected values of the spin must be calculated for whole lattice sites in every MC step, more CPU time seems to be wasted comparing to the Ising model, this shortcoming can be overcomed by technical means. Analysis Hamiltonian (4) found that the eigen-state of spin S i is only dependent on the field )space before the simulation, it is called directly in the later calculation. In this way, it doesn't take more time than simulation of Ising model with same size. In addition, in order to improve the calculation accuracy and efficiency to solve the Hamiltonian (4), we can transform the Hermitian matrix into a real symmetric matrix by rotation of coordinate system, which simplifies the calculation and improves the calculation accuracy.
Using the selected state, we can calculate the thermodynamic quantities of system by the following formule. Defining the a operator,  periodic boundary condition is used. For a definite temperature, 3×10 4 MC steps per spin are performed to equilibrate the system, succeeded by a run of 5×10 4 MC steps per spin for average the physical quantities. For brevity, we take the Boltzmann constant k B as 1. Figures 1 and 2 show the temperature dependence of the thermodynamic properties of ferromagnets with SC, BCC and FCC structures for S=1/2 and 1 ferromagnets respectively. With temperature decreases from high temperature region to a definite point, the magnetization increases suddenly, an inflection point appears on energy versus temperature curve, and the specific heat and the susceptibilities show peaks at this point. It is obvious that phase transition from paramagnetic phase to ferromagnetic phase occurs in the system and temperature point at which specific heat and susceptibility get peaks is defined as Curie temperature or critical temperature. It is also found that the decrease rate of specific heat on both sides of the transition point is different, above the critical point, the specific heat decreases faster with the temperature, showing the characteristics of the second order phase transition. Although we apply a small magnetic field to the system in simulation (h/J=0.001), the positions of the peak points for specific heat and susceptibilities are almost the same. There are obvious fluctuations in the magnetizations for S=1/2 ferromagneties below the critical point ( figure 1(a)), which leads to the fluctuations of longitudinal and transverse susceptibilities at the left side of the critical point. Above the Curie temperature, the values of longitudinal and transverse susceptibilities are almost same for D=0, which is the natural result of the MMFMC. When D>0, the longitudinal correlation of spins is stronger, the longitudinal susceptibility is bigger than the transverse susceptibility at the same temperature ( figure 2(e)). At zero temperature, the thermodynamic quantities calculated are the exact results of the ferromagnets. The magnetization is their classical saturation value, the ground state energy shows the parallel arrangement of spins, and the specific heat and the susceptibilities take zero value. In order to verify the validity of the calculation results under the present size, we changed the sizes of the simulation cell for S=1 system, and found that the thermodynamic properties and Curie temperatures of SC, BCC, and FCC systems changes little as  L L L a ,

Numerical results and discussions
( ) which shows that the current simulation size in present method is suitable for simulation of the 3D Heisenberg system. Table 1 give the Curie temperatures of Heisenberg ferromagnets with SC, BCC and FCC structures for several quantum number S. For comparison the results obtained by THE, QMC and GFM are also given in the table. It is found that the Curie temperatures obtained by present method is consistent with those obtained by the other methods. The Curie temperature, Tc/J is 0.836, 1.227, and 1.964 for S=1/2 ferromagnets with SC, BCC and FCC structures, respectively, and they are 0.839, 1.260, 2.008 in THE [19,20]. QMC simulation gives 0.839 [24,26] for SC ferromagnet. The present results are slightly smaller than those of THE and QMC. The results of the pseudofermion functional renormalization group (PFFRG) method for BCC ferromagnet with S=1/2, 1.45(one-loop) and 1.37(two-loop) [27], also bigger than the present result. For S=1 ferromagnet with three structures, the Curie temperatures obtained with the MMFMC method are 2.701, 3.901, and 6.180, respectively, and they are 2.599, 3.788 for SC and BCC structures respectively in THE [20]. For S>1, the Curie temperatures obtained by the present method are little higher than those obtained by THE and GFM [10,11,16,20], but they are very close to the THE results for BCC lattice. Tc/[JS(S+1)] increases slowly with spin quantum number, which is agreement with that obtained by THE. In GFM with RPA approximation (10,11), the Curie temperature can be expressed as where -F 1 ( ) is a structure only factor, so the Curie temperature is independent of the spin quantum number. The introduction of the higher-order Green functions may improve the above results and the Curie temperature was found to increase with the spin quantum number [16,28], and the specific heat can be also calculated for the low-dimensional Heisenberg systems [29]. From table 1 we also found that the present MMFMC method can be applied to any spin system compared with the QMC method, and the amount of calculation is not very large.
We give the Curie temperature as a function of anisotropic constant in figure 3. It is found that the Curie temperature increases slowly with the anisotropic constant, which is consistent with the dada given in [15] and [16]. In the transition from paramagnetic to ferromagnetic phase, the quantum state with the extreme spin value is the most possible due to the ferromagnetic interaction between spins. With the introduction of the easy-axis anisotropy, the quantum state will be further strengthened, but the essence remains unchanged, so the Curie temperature of system increases slowly with the anisotropic value. The crossover of the Curie temperature versus anisotropy curves for different spin S is due to the calculation errors. In MMFMC method, only the eigen equation of a single spin in a field needs to be solved, so the existence of the single-ion anisotropic term does not bring about more computation.
In MFT for ferromagnet with S=1, the Curie temperature satisfy the following relation, It can be found that the Tc increases smoothly with the anisotropic value D. Although the quantum characteristics of a spin in a micro mean-field is considered, the correlation of the nearest-neighbor spins is decoupled, the quantum character between spins is omited, so we cannot get the different critical temperatures for ferromagnet and antiferromagnet with the same structure and spins in MMFMC like the other methods [20,25,30].