Spin waves in one dimensional magnetic material

Using Heisenberg model, the equations of motion for the dynamic properties of spin waves in three dimensions were obtained and solved analytically up to an exponential operator representation. Second order Suzuki Trotter decomposition method with evolution operator solution was applied to obtain the numerical solutions by making it closer to real spin systems. Computer based simulations on systems in micro canonical ensembles in constant-energy states were used to check the applicability of this model for one dimensional lattice by investigating the occurrence, temperature dependence and spin-spin interaction dependence of the spin waves. A visualization technique was used to show the existence of many spin wave components below the Curie temperature of the system. In the magnon dispersion curves all or most of the spin wave components could be recognized as peaks in the dynamic structure factor. Energy conservation of the algorithm is also shown.


INTRODUCTION
Spin waves are propagating disturbances in ordering of spin orientations of magnetic materials. Due to spin-spin interactions between nearby particles, spin which is an intrinsic property of a quantum mechanical particle is a dynamic property. Understanding spin dynamics is important because they have many applications in nanotechnology and quantum computing. Existence of spin waves can be experimentally observed by using various inelastic scattering methods such as scattering of neutrons, light, electrons and ferromagnetic resonance. However, to understand their occurrence, dynamical equations corresponding to quantum spins have to be solved. A realistic model would be an array of atoms on a lattice. The major difficulty in solving such a system analytically is that due to the existence of spinspin interactions it is usually chaotic and too complex.
Even though the spin is entirely a quantum mechanical concept, classical models that describe the dynamics of spin systems exist. Using these models, spin systems can be studied analytically or numerically. A realistic classical spin model should represent a spin as a proper vector in three dimensions. Though most classical models like Ising model can be analytically solved, they do not represent spin as three dimensional vectors. Classical Heisenberg model can be used to study the spin dynamics in three dimensional vector representations. But for a three dimensional lattice, the Heisenberg model cannot be solved analytically. Landau et.al [1] studied magnetic excitations at critical and low-temperature using simple classical spin models with true dynamics, governed by equations of motion. One major limitation in spin dynamical simulations is the spatial and temporal constraints due to the finiteness of the computational resources. Although to speed up the simulation it is desirable to choose the  LTD, 2015 largest possible time step size, according to Landau et.al, the size of the time step is severely limited by the accuracy within which the conservation laws of the dynamics are obeyed. Total energy and net magnetization are good examples for conserved entities according to the isotropic Heisenberg model. They proposed the application of the Suzuki-Trotter method as an alternative and the single most important quantity that is extracted from the numerical results is the dynamical structure factor. In order to compute the dynamical structure factor, however, space and time displaced spin-spin correlations must be first computed. In order to compute the spin-spin correlations, Tao et.al [2] described an estimating procedure for the ensemble averages using a large number of starting configurations and used a technique that avoided storage of a huge amount of data to verify the existence of spin waves below Curie temperature. In most of the studies one level of nearest neighbours were coupled using a single interaction parameter or dimensions were limited to two. Costa et.al [3] using Monte Carlo methods and molecular-dynamic techniques simulated two-dimensional Heisenberg model with anisotropy. Due to added complexity of Heisenberg model that makes analytically solving the system extremely difficult if not impossible for some cases, the study of spin waves is mostly based on Ising model and XY model.
In the present work, a model based on Heisenberg model was developed to observe the dynamic property of the formation of spin waves in magnetic material. The simulations were carried out up to the second nearest interaction parameters. Allowing energy transfer with the environment is difficult in a simulation because it makes the system extremely complex. Therefore it is convenient to simulate the system in a micro canonical ensemble where the system is in a constant-energy state. A more realistic system would be in constant temperature rather than constant energy. To initiate the dynamic simulation, it was first necessary to bring the spin system into thermal equilibrium with a heat bath of known temperature. This requirement was achieved by using Metropolis algorithm, a mainstay of Monte Carlo methods and statistical simulations. In the program, a set of simulation units were used rather that SI units where the Boltzmann constant was set to unity; the interaction parameter was absorbed into simulation time; the second interaction parameter was defined relatively to the first interaction parameter; the gyromagnetic factor was taken as unity and the temperature was measured relative to Curie temperature of each spin system.

THE DYNAMICS OF THE SPIN SYSTEM
The time evolution of a system that can be modelled as being exactly in one state at a given time is described by the master equation. The probability that the system will be in state at time t is () t   and the rate of change is given by the master equation as, R   is the transition rate from state to  and the total probability of the system for all the states is one. The expectation value of any thermodynamic variable Q at time t can be calculated as (2) with the transition probabilities satisfying sum rule, For a canonical ensemble, the equilibrium probabilities ρ μ at temperature T are given by, where E μ is the energy associated with the state μ, k B is the Boltzmann constant and the partition function Z is given by, Therefore, at equilibrium the expectation value of the thermodynamic variable is This equation can be rewritten in terms of selection probability   g   , the probability that the algorithm will generate a new target state ν from an initial state μ and acceptance ratio   A   , the probability that the selected state will be accepted, In the simulation, acceptance ratio should be made as close to unity as possible. Assuming that there are M discrete states in the phase space of the spin system such that the system can be in each of the possible states with the same probability that is with uniform probability and thereby setting  

ILCPA Volume 47
The largest of the acceptance ratios is made to be equal to one while the others are adjusted to satisfy the constraint A fully quantum mechanical treatment of spin dynamics is extremely difficult and inefficient for numerical simulations. Therefore, the spin system is treated as a classical system of interacting vectors on a lattice governed by a three-vector model, known as classical Heisenberg model. The lattice is considered as a three-dimensional grid of lattice sites. Since the magnitude of spin for each lattice site is equal, it is convenient to represent the classical spin corresponding to lattice site i by a unit vector 1 i 

S
. Classical isotropic Heisenberg model can be described by the Hamiltonian, where J is the spin-spin interaction parameter and i S is a classical spin located on a lattice site i. Dynamics of the Heisenberg model is governed by the equation of motion [4]: where γ is the gyromagnetic factor and eff i H is the effective field given by, NN(i) denotes the set of nearest neighbouring lattice sites of i site. Equation 12 is really a set of N equations for 1, 2,3..., . iN  Because the effective field has units of interaction parameter J, the time variable used in the spin dynamics simulations is also measured in units of J [5] .
If the classical spin and effective field vectors are denoted by matrices, then by setting 1   , equation 12 can be rewritten in matrix form Assuming continuous time and R is a constant, the solution to the above differential equation takes the form  ˆR t y t Ce  . However even if time is discrete and R is a matrix for discrete time steps τ, a solution can be written as [6], In the above equation, the matrix operator . For any skew-symmetric matrix R, R e is a rotation matrix [7]. Therefore the evolution operator in equation 16 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 15 can be solved using Suzuki Trotter decomposition method [8]. When A and B are arbitrary non commuting operators, the product of the exponential operators is regarded as an approximate decomposition of the exponential operator with correction terms of the second order of τ. A second order symmetric Suzuki Trotter decomposition is, Considering a lattice space with spin vectors corresponding to each lattice site, the space-displaced time-displaced spin correlation function can be defined as [1]: is the k-component of spin vector S located at site r at time t in the lattice space. The closed 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 ( , ) S  q and it is related to the static structure factor () S q by the following integral equation [9]: where q is the scattering vector or wave vector and ω denotes the frequency. This can be computed by a Fourier transformation of the spin-correlation function: The relation of the static structure factor and dynamic structure factor given in equation 22 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 23 can be rewritten as,

SIMULATION OF SPIN DYNAMICS
The spin system was simulated using MatLab 7.14 syntax. The Metropolis algorithm was used to obtain thermal equilibrium and the Suzuki Trotter decomposition was used to solve the dynamics of the system. To achieve the thermal equilibrium by heat exchanges with a constant temperature environment, the lattice was initially defined with each spin oriented in z direction. Then for each lattice site, the spin orientation was randomized by using spherical polar coordinate system in order to avoid the normalizing of spins each time a randomization was done. The spin randomization was achieved by generating random numbers for the for ϕ. However for an extremely large set of Monte Carlo iterations, the spin system was not in proper equilibrium, as this method does not pick random points on the surface of a unit sphere properly ( figure 2(a)). Anisotropy in z direction is clearly visible in the plot. In order to avoid the anisotropy, the Marsaglia's method [9] was adopted to generate random orientation in spin in Monte Carlo simulation. Figure 2(b) shows 1000 randomly generated spin orientations using Marsaglia's method. Using the new energy term, the Metropolis algorithm was applied.
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 is done in micro canonical ensemble. To get the average over canonical ensemble, about 50 simulations per each system were performed and their correlation functions were averaged out. A file with two variables including the correlation function matrix and the simulation count was created. Each time a new simulation was performed, correlation data were added to the matrix and the simulation count was increased by one. At a later stage these data were used to compute the averages. After computation of the integrand of the dynamic structure factor obtained in equation 24, the numerical integration was done using trapezoidal rule. By setting magnetization direction as the +z direction, so as to identify the spin waves, the transverse components in

ILCPA Volume 47
A method was developed to implement the spin dynamic simulations, with second nearest neighbour spin interactions on one dimensional lattice. The first nearest neighbours of i th spin are coupled to it via the usual interaction parameter, whereas second nearest spins are coupled to the same spin with a much smaller secondary interaction parameter as they are shielded by first nearest neighbours.

SPIN WAVES IN ONE DIMENSIONAL LATTICES
In one dimensional lattices, only one level of nearest neighbours were considered assuming that all the neighbours of a spin S i interact with it via the same interaction parameter. A lattice space with N sites considered is shown in figure 3(a). In order to show the lattice decomposition of the computational algorithm clearly, the atoms are denoted as cubes in nonstandard notations. In the simulation, the nearest neighbours of site 2 were chosen to be site 1 and 3 etc. Since only the nearest neighbours were taken into account when the effective fields were computed, for the effective fields H effAi to be essentially constant while a rotation was performed in the sub lattice A, the nearest neighbours for site 1 and site N, circular boundaries were used so that site 1 has site 2 and site N as its neighbours, while site N has site N-1 and site 1 as its neighbours. If this decomposition is not used, the effective field does not stay fixed while rotating the spins since all the spins must be concurrently updated. Suzuki Trotter decomposition together with above sub lattice choice makes this rotation possible.
For the one dimensional spin chain, the nearest and second nearest neighbours of i th spin i S are shown in figure 4. To integrate such a spin system, two sub lattices are inadequate. When two sub lattices are used, as the evolution operator which depends on the second nearest neighbours cannot be applied on i th spin within discrete time step τ, the evolution operator itself changes as it is on the same sub lattice that i th spin lies. To eliminate this problem, the sub lattice decomposition shown in figure 5 was proposed. One dimensional, spin systems were simulated for systems with single coupling constant and two coupling constants separately and for different temperatures. Table 1 gives the data for the simulation. At least 50 simulation runs were made and for most systems, a single run produces few megabytes of data and therefore listing all the data here is impossible.

ILCPA Volume 47
At 0.1 C T (T C is the Curie temperature), a system of 30 spins was simulated for nearest neighbour coupling. The Monte-Carlo update of the system is shown in figure 6. Metropolis algorithm was used to obtain the thermal fluctuations. The plot ensures that the system is in required temperature. Energy per spin has stabilized around 0.9  despite some stochastic fluctuations. Once the integration is done, stored data can be used to identify spin waves. To visualize the spin orientations in the direction normal to magnetization, a technique was developed. The x component of spin orientations plotted against time is shown in figure 7. This graph is possible only for one-dimensional systems and shows very interesting information about the time development of the system. The figure shows clearly the wave-like propagations present in the spin system. Diagonal patches of similar colours suggest that spin orientations tend to propagate like a wave as time develops. Furthermore, two sets of diagonal lines are clearly visible. One set correspond to larger patches and the other set shows smaller and denser lines. Those smaller lines are steeper and more densely packed than the larger patches suggesting they correspond to faster propagating waves and smaller wavelengths. That is they correspond to spin wave components of larger propagation vectors. Larger patches occupy larger number of spins across the lattice. Since they are not as steeper as smaller lines, they propagate slowly. The propagation of spin waves is illustrated in the figure 8. Spin orientations change with respect to time influenced by nearby spins, giving rise to spin waves. The noticeable blue-coloured disturbance (left) has propagated leftwards as time develops.

STRUCTURE FACTOR
The computation of structure factor gives more components of waves that are not quite visible in the figure 7.

ILCPA Volume 47
The computed dynamic structure factors against frequency The dynamic structure factor obtained for different propagation vectors confirmed the existence of different components of spin waves corresponding to different momentum transfers. Though for dynamic structure factor, smooth curves should be expected, many oscillations are present in all graphs due to the extremely small number of atoms in the system. Another reason for these oscillations is the computation of ensemble average. In reality, even 30 atoms system can occupy infinitely many equilibrium states in the canonical ensemble. Average over this infinitely many ensembles is approximately computed by only using 50 or 100 simulations. Thus the smaller peaks resemble the noise in the data. To compare the peaks, it is useful to plot the Magnon dispersion curve (figure 10) which shows all the possible peaks in a single graph. A series of spin waves corresponding to different energy transfers and different momentum transfers is observable in the figure. The graph shows that for the spin waves corresponding to higher momentum transfer, the frequency gaps become smaller and smaller.
After the system was brought to thermal equilibrium, the dynamical equations of the system were used to simulate the time development. If the integration method was accurate enough energy of the system or energy per spin should stay constant at the final energy state after the Monte Carlo simulation. This can be checked by plotting the energy versus simulation step presented in figure 11. The energy is conserved in a one-dimensional system with a single interaction parameter. In the integration phase, it is clear that Suzuki Trotter algorithm very sharply preserving the energy of the system. This is indicated by the horizontal line corresponding to energy per spin during the integration steps. Another measure is to check whether the magnetization of the system is preserved by the integration method. Thus the magnetization of the spin configuration was computed and stored in each integration step. Figure 12 shows the magnetization of the one-dimensional system against time. Many oscillations in the graph shows that second order Suzuki Trotter algorithm does not preserve magnetization of the system exactly even for the one-dimensional lattice. However, the variation of magnetization is extremely small (10 -4 ). The magnetization stays within +0.7674 and +0.7675 during the simulation. Since in the used temperatures all spin wave components show much larger oscillations, the effect of these small fluctuations can be neglected in this study of one-dimensional systems.

ILCPA Volume 47
Random spin orientations can be expected if a plot is generated for a system above its Curie temperature T c . Figure 13 shows the time development of the one-dimensional system at a temperature above its T c . At this temperature, there is no evidence of patterns of lines traced by the spin waves. Even though the graph looks completely random it is not random, but a result of an extremely complex and chaotic trajectory traced in a constant energy surface of the micro canonical ensemble. Figure 14 shows the path of a point on the unit surface onto which a single spin of the system points as time develops.

International Letters of Chemistry, Physics and Astronomy Vol. 47
One-dimensional spin system was simulated at two different temperatures 0.1T c and 0.2T c to find out the temperature dependence of the spin waves at these temperatures. These are shown in figure 15. They show that the dynamic structure factors corresponding to the same momentum transfer at both temperature. Two peaks correspond to the two different temperatures which are both below the Curie temperature T c of the system. According to the figure, the spin waves corresponding to a momentum transfer of 11π/15 tends to transfer less energy as temperature is increased. At 0.1 T c the spin wave corresponds to a frequency of 3.2325 rad.s -1 whereas at 0.2 T c the frequency decreases to 3.0079 rad.s -1 . That is, the corresponding wave component decreases in frequency as temperature is increased. The peak corresponding to higher temperature looks noisier in the plot. However both systems were simulated hundred times and their results were averaged out. The temperature dependence holds for other momentums as well. The red peak is shifted leftwards relative to corresponding black peak which are at a lower temperature. Therefore spin waves corresponding to a given propagation vector tend to transfer less energy as temperature increases. That is, magnons with the same momentum tend to decrease their energies as temperature is increased. Also the frequency shifts Δω becomes larger with momentum transfer.

CONCLUSIONS
The second order Suzuki Trotter decomposition technique along with evolution operator solution to spin equations of motion is successful in numerically solving one-dimensional spin systems and can be successfully extended to include a second nearest interaction parameter. Classical Heisenberg model produces extremely complex and chaotic motion for spin orientations. A visualization method was used to show the existence of many spin wave components below the Curie temperature of the system. Even though two clearly spin wave components were visible, a careful observation shows many other spin wave components corresponding to different propagation vectors and energies. To show the existence of these less obvious components, the correlation functions were computed and their space-time Fourier transformations were plotted. In the magnon dispersion curves all or most of the spin wave components could be recognized as peaks in the dynamic structure factor. The spin waves occur for spin systems below Curie temperature in a group and various components of these correspond to different momentums and energies. As temperature is increased but kept below Curie temperature, the frequency of the spin wave components decreases regardless of the propagation vector. There is a strong relation between spin-spin interaction and energy transfer of the spin waves. With the availability of more interactions, the spin waves can transfer more and more energy. Energy was conserved in the algorithm.