Spin waves in two and three dimensional magnetic materials

The equations of motion for the dynamic properties of spin waves in three dimensions were obtained using Heisenberg model and solved for two and three dimensional lattices analytically up to an exponential operator representation. The second order Suzuki Trotter decomposition method was extended to incorporate second nearest interaction parameters into the numerical solution. Computer based simulations on systems in micro canonical ensembles in constant-energy states were used to check the applicability of this model for two dimensional lattice as well as three dimensional simple cubic and bcc lattices. In the magnon dispersion curves all or most of the spin wave components could be recognized as peaks in the dynamic structure factor presenting the variation of energy transfer with respect to momentum transfer of spin waves. Second order Suzuki Trotter algorithm used conserved the energy.


INTRODUCTION
Even though the spin is entirely a quantum mechanical concept, using classical models spin systems can be studied analytically or numerically. Most classical models like ising model do not represent spin as three dimensional vectors. But the classical Heisenberg model uses the three dimensional vector representations for spins. Spin-dynamics technique was used by Landau et.al [1] to study the critical and low-temperature magnetic excitations using simple classical dynamical spin model governed by equations of motion. According to the isotropic Heisenberg model, the total energy and net magnetization are good examples for conserved entities. To investigate the existence of spin waves in paramagnetic bcc iron, Tao et.al [2] used computer simulations based on classical Heisenberg model on a cubic lattice with periodic boundary conditions and obtained the equilibrium states by applying Monte Carlo method. In most of the studies one level of nearest neighbours were coupled using a single interaction parameter or dimensions were limited to two. To simulate a twodimensional Heisenberg model with anisotropy, Costa et.al [3] used Monte Carlo methods and molecular-dynamic techniques.
In the present work, a model based on Heisenberg model used by us to study the spin dynamics of one dimensional magnetic material [4] was extended to observe the dynamic property of the formation of spin waves in two dimensional lattices as well as three dimensional simple cubic and body centred cubic lattices of magnetic materials. The simulations were carried out up to second nearest interaction parameters. The systems were considered to be in micro canonical ensembles in constant-energy states. To initiate the dynamic simulations, the spin systems were brought into thermal equilibrium with heat baths of known temperature by using Metropolis algorithm which is a mainstay of Monte Carlo methods and statistical simulations. In the program a set of simulation units were used rather than SI units where the Boltzmann constant was set to unity, interaction parameter was absorbed into simulation time, second interaction parameter was defined relatively to first interaction parameter, gyromagnetic factor was taken as unity and temperature was measured relative to Curie temperature of each spin system.

THE DYNAMICS OF THE SPIN SYSTEM
The time evolution of the probability () t   of a system in state at time t is given by the master equation: where () R   represents the transition rate from state  to . The sum of the probabilities of all the states is one. For a canonical ensemble, the equilibrium probabilities ρ μ ( ( ) ) t    at temperature T are given by, where E μ is the energy associated with the state μ, k B is the Boltzmann constant and Z is the partition function. Therefore the expectation value of the thermodynamic variable Q is The condition for thermal equilibrium 0 d dt    leads to an equation that can be rewritten in terms of selection probability g(μ→ν) and acceptance ratio A(μ→ν). Assuming that there are M discrete states in the phase space of the spin system and that the system can be in each of the possible states with the same probability (i.e. a uniform probability) and thereby setting g(μ→ν)=1/M, The largest of the acceptance ratios is made to be equal to one while the others are adjusted to satisfy the constraint The exponentiation of the matrix operator R can be identified as another matrix operator given by . For any skew-symmetric matrix R, R e is a rotation matrix [7]. Therefore the evolution operator in equation 8 rotates each spin in () yt about the effective fields by an amount proportional to the field magnitude. The set of ordinary differential equations derived in equation 8 can be solved using Suzuki Trotter decomposition method [8]. If the lattice is divided into four interpenetrating sub lattices denoted by A, B C and D, then the effective Hamiltonian The space-displaced, time-displaced spin correlation function and its space-time Fourier transform were used to study of spin dynamics. Considering a lattice space with spin vectors corresponding to each lattice site, the space-displaced time-displaced spin correlation function is defined as [1]: is the k-component of spin vector S located at site r at time t in the lattice space. The angle brackets 〈〉 denote the average over the ensemble. The information about inter-particle correlations and time evolution of the system is contained in the dynamic structure factor , () The static structure factor () k S q can be obtained by integrating the dynamic structure factor over all the frequencies [9]. Dynamic structure factor can be computed by applying space-time Fourier transformation of the space-displaced, timedisplaced spin-correlation function: The relation of static structure factor and dynamic structure factor given in equation 14 shows that the dynamic structure factor have to be a time translation invariant and space translation invariant quantity. Therefore by taking only the real parts of the exponentials, the equation 14 can be rewritten as,

THE SPIN DYNAMICS SIMULATION
For an original lattice with each spin oriented in z direction, the Metropolis algorithm was used to obtain thermal equilibrium by heat exchange with a constant temperature environment. In achieving thermal equilibrium with a constant temperature environment, for each lattice site, the spin was randomized using spherical polar coordinate system. To the new energy term, the Metropolis algorithm was applied. In order to avoid the anisotropy, the Marsaglia's method [9] was adopted to generate random orientation in spin in Monte Carlo simulation. This is shown in figure 2. The Suzuki Trotter decomposition was used to solve dynamics of the system. The space-displaced, time-displaced spin-correlation function ,) ( k C ' t  rr was computed and stored as a function of '    r r r and time t. For a single simulation run however, the spin system is in constant energy since the simulation was done in micro canonical ensemble. To get the average over canonical ensemble, about 50 simulations per each system were performed. A correlation function matrix and the simulation count were created. After the computation of the integrand of the dynamic structure factor in equation 15, the trapezoidal rule was used for the numerical integration. By setting magnetization direction as the positive z direction, the transverse components in x and y directions were identified as () , ) , ( xy SS   qq .

SPIN DYNAMICS OF TWO DIMENSIONAL LATTICE
A method was developed to implement spin dynamic simulations with second nearest neighbours on two dimensional and three dimensional simple cubic lattices. The first nearest interactions of i th spin were coupled to it via the usual interaction parameter, whereas second nearest spin interactions were coupled to the same spin with a much smaller secondary interaction parameter. A two dimensional lattice space with N L L  sites was constructed as shown in figure 3. The atoms are denoted as cubes in non-standard notations to show the lattice decomposition of the computational algorithm clearly. For a two dimensional lattice, there are four nearest neighbour sites per each site. In the simulation, neighbours of second site were chosen as the first and third site etc. Circular boundaries were used for the two end sites, so that site 1 possesses site 2 and site N as its neighbours, while site N possess site 1 N  and site 1.
The decomposition of the sub lattice was done by assigning each odd-numbered site to sub lattice A and each even-numbered site to sub lattice B. The resulting lattice definition is shown in figure 4. Each spin in sub lattice A (red), experience the effective fields due to its nearest neighbours which are all formed from sub lattice B (orange). Therefore when performing a rotation on sub lattice A, the effective fields due to sub lattice B stay fixed and vice versa. On a two dimensional lattice space, the way the first and second nearest neighbours were defined is shown in figure 5. To integrate such a spin system, two sub lattices are inadequate. If only two sub lattices are used, evolution operator cannot be applied on i th spin since within discrete time step , the evolution operator itself changes since it depends on secondnearest neighbours which are on the same sub lattice that i th spin is in. To eliminate this problem, the sub lattice decomposition shown in figure 6 was introduced. As a result of this sub lattice definition, all the first-nearest neighbours and second-nearest neighbours of any spin S i are due to other sub lattices and thereby gave a successful decomposition.
To simulate the two dimensional spin systems the data were collected for systems with single coupling constant and two coupling constants separately and for different temperatures and at least 50 simulation runs were made. These data are tabulated in table 1. For most systems, a single run produces few megabytes of data and all the data is not listed. A two-dimensional lattice with 100 sites 10 10  was simulated at a temperature of 0.5 C T . Monte Carlo updates simulating equilibrium and energy fluctuations are shown figure 7. The Monte Carlo simulation shows a more relatively stable system than for the one dimensional case for the same temperature. In the one dimensional system simulation used only 30 atoms whereas for two-dimensional system 100 atoms were used [4]. With the increase in number of atoms in the simulation, the system's thermal equilibrium becomes much stable and energy per spin shows fewer fluctuations. To find out the effect of nearest neighbour coupling on energy per spin at thermal equilibrium, the same system was simulated considering first and second neighbour interactions. The resulting energy fluctuation curve is shown in figure 8. For only the first nearest interactions, the equilibrium state of energy per spin is around 1.5  . At the same temperature, when both first and second nearest spin interactions were included in coupling, the system is in a much lower equilibrium state of energy per spin ( 2.1)  . The dynamic structure factors can be computed and compared of systems at the same temperature up to the first and second nearest neighbour spin interactions to see the dependence of nearest neighbour coupling and the spin wave frequencies. Figure  9 shows the resulting dynamic structure factors of twodimensional systems with same lattice size simulated at the same temperature and for the same wave vector. For the same propagation vector, the second neighbour coupled system transfers considerably greater amount of energy than the nearest neighbour coupled system. In this particular case, second-rank coupled system transfers almost twice the energy of what first-rank system does. The system used to generate above graphs was a two dimensional array of 10 10  spins. Since under these conditions only five momentum transfers are possible, they were plotted in figure 10 to find out if this pattern continues. These figures establish the expected dependence of energy transfers and spin interactions. With the availability of more interactions, the spin system tend to transfer more energy at the same temperature relative to T C , suggesting that when more spins are available to interact, ferromagnetic spin systems generate spin waves of higher frequencies.

SPIN DYNAMICS OF THREE DIMENSIONAL SIMPLE CUBIC LATTICE
In the case of the three dimensional lattices, the spin simulation were performed for a simple cubic (sc) crystal structure and a bodycentred cubic (bcc) crystal structure up to second nearest spin interactions. For three-dimensional lattice with simple cubic crystal structure, numbering lattice sites and decomposing into sub lattices is straightforward and similar to the previously discussed twodimensional situations. In three dimensions, there are six nearest neighbours per spin. For a spin in a given sub lattice, all the neighbours come from the other sub lattice (Figure 11). To integrate the simple cubic spin system, the decomposition into eight sub lattices as shown in figure 12 was proposed. For a spin i S belonging to a particular sub lattice, all the first and second-rank nearest neighbours come from other sub lattices. Therefore evolution operator can be decomposed to apply as a series of operators on spin i S . The systems of which spin wave occurrence is investigated are tabulated in table 2 and for each system at least 50 simulation runs were made. Using MatLab a 666  system was simulated successfully. The spin energy variation in 500 Monte-Carlo steps for a simple cubic lattice with only nearest neighbour interactions are shown in figure 13. A three dimensional system usually having more spins, 216 spins in the above particular case arrive at the stable equilibrium with even less fluctuations than two dimensional lattices. The magnon dispersion curve showing all the spin components can be generated for this system by computing the dynamic structure factor for the individual momentum vectors and plotting them in the same graph. This is shown in figure 14. The magnon dispersion curve with first-rank coupling suggests that energy gaps between magnons close up as their momentums increase. In the above system the peaks correspond to 1 0.911 rad.s ,  1 2.768 rad.s  and 1 3.683 rad.s  . Therefore, spin waves corresponding to consecutive momentums show smaller energy differences at higher momentums.  Figure 12: Sub lattice definition of sc lattice with two coupling constants Figure 13: Thermal equilibrium of a sc lattice system with one coupling constant at 0.5 T C A three-dimensional simple cubic (sc) lattice with a lattice decomposition shown in figure 12 was simulated with the second coupling constant 0.3 times that of the first coupling constant. The spin energy of the sc system with these two interaction parameters for a 0.0001 Suzuki Trotter integration step size is shown in figure 15. According to the figure, the energy conservation can be maintained for a sc spin system. However it requires extremely small step sizes. For such extreme step sizes, generating data enough to visualize or compute space-displaced, timedisplaced spin correlation functions is therefore extremely tedious.

SPIN DYNAMICS OF BODY CENTERED CUBIC LATTICE
For a body centered cubic structure (bcc), the first and second nearest neighbour definition is slightly complex. First, a numbering method should be applied to identify the spins. In doing so, the lattice is defined as two interpenetrating simple cubic like lattices as shown in figure 16. For one lattice, the usual numbering system is applied while for the other lattice, a similar numbering system is applied but starting from the largest number in first lattice as 1  . From the structure in the figure 16, two ranks of nearest neighbours are identified. These are shown in figure 17. For a bcc structure, the lattice decomposition is much easier than for a simple cubic structure. It can be done by just decomposing the two simple cubic like lattices each into two sub lattices where a sc lattice was decomposed with nearest neighbours as shown in figure 18 with a single coupling constant. With proper lattice site numbering, sub lattice of any given spin can be identified by its index being odd or even just as in the two dimensional lattices. For the bcc structure, therefore, only four sub lattices result in and there are fourteen nearest neighbours including eight first nearest neighbours and six second nearest neighbours. The systems of which spin wave occurrence is investigated are tabulated in table 3. For each system at least 50 simulation runs were made. A bcc spin system defined on a 10    The Monte Carlo simulation yields a much stable equilibrium. According to the figure, Suzuki Trotter integration preserves energy almost exactly as suggested by the flat curve to the right side of the plot. Therefore above figure suggests that decomposition method works successfully.

CONCLUSIONS
This study was focused at utilizing computational power to study complex dynamics of magnetic systems, mainly studying the propagating disturbances in the ordering of spin orientations. Second order Suzuki Trotter decomposition technique along with evolution operator solution to spin equations of motion was successful in numerically solving two-dimensional and three-dimensional spin systems and can be successfully extended to include a second interaction parameter for two and three dimensional systems with simple cubic and body centred cubic structures. Using the developed MatLab code, many systems were simulated and analysed to identify spin waves and their properties in various conditions. Energy conservation of the algorithm was shown for each situation. For a simple cubic lattice with two interaction parameters, there was a significant energy drift. However, this can be attributed to the integration method used. Second order Suzuki Trotter algorithm is accurate only to the second order of step size. However the algorithm used was conserving energy when extremely small integration step sizes were used for very short time periods. In the magnon dispersion curves all or most of the spin wave components could be recognized as peaks in the dynamic structure factor. They clearly presented the variation of energy transfer with respect to momentum transfer of spin waves.
Classical Heisenberg model produces extremely complex and chaotic motion for spin orientations. Spin waves occur in one-dimensional, two-dimensional, and three-dimensional spin systems below Curie temperature, and when they occur, they occur in a group and various components of them correspond to different momentums and energies. As temperature is increased but kept below Curie temperature, frequency of spin wave components is decreased regardless of the propagation vector. There is a strong relation between spin-spin interaction and energy transfer of spin waves. With more interactions available spin waves can transfer more and more energy.