Theoretical analysis of neutron scattering results for quasitwo dimensional ferromagnets
Abstract
A theoretical study has been carried out to analyze the available results from the inelastic neutron scattering experiment performed on a quasitwo dimensional spin ferromagnetic material . Our formalism is based on a conventional semiclassical like treatment involving a model of an ideal gas of vortices/antivortices corresponding to an anisotropic XY Heisenberg ferromagnet on a square lattice. The results for dynamical structure functions for our model corresponding to spin, show occurrence of negative values in a large range of energy transfer even encompassing the experimental range, when convoluted with a realistic spectral window function. This result indicates failure of the conventional theoretical framework to be applicable to the experimental situation corresponding to low spin systems. A full quantum formalism seems essential for treating such systems.
PACS: 78.70.Nx–inelastic neutron scattering in condensed matter, 75.10.Jm–Heisenberg model,
75.30.Kz–KosterlitzThouless transition in magnetic systems
1 Introduction
Low dimensional and in particular two dimensional magnetism has attracted a great deal of interest in the past three decades [1, 2, 3]. In particular, in one dimension the existence of both solitonic and spin wave excitations were thoroughly studied through inelastic neutron scattering experiments as well as theoretical analysis for [4]. Similar studies were carried out searching for topological excitations in various quasione dimensional systems which are almost ideal realization of nearest neighbour Heisenberg antiferromagnetic chain [5].
In many of the above systems the experiments showed central peak (peak corresponding to ) in the dynamical structure function when plotted in constant scan. This motivated the experimentalists further to investigate two dimensional and quasitwo dimensional magnetic materials. With the availability of improved quasi twodimensional ferromagnetic and antiferromagnetic materials investigations along this line has become todays one of the primary interests both theoretical and experimental. These include layered systems such as , , magnetically intercalated graphites such as , layered ruthenates, layered manganites and high cuprates [1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12]. Moreover large amount of information on the spin dynamics, extracted from inelastic neutron scattering are available. Advances in numerical and computational techniques have also contributed to the understanding of both spin wave and topological excitations [13, 14, 15, 16, 17].
On two dimensional magnetic systems the concept of topological order was introduced by Kosterlitz and Thouless and independently by Berezinskii [18, 19]. Their ideas backed by analytical and numerical calculations led to the proposal for the existence of topological vortices and antivortices in a typical ferromagnetic XY model on a two dimensional lattice. According to these ideas vortices and antivortices are frozen as bound pairs below certain transition temperature called or and above this temperature they become mobile and nearly free [18, 19].
In this work we initiate a theoretical investigation regarding the applicability of a semiclassical like treatment of the dynamics of topological excitations to the inelastic neutron scattering results for real systems [6, 13, 14, 20]. In recent years inelastic neutron scattering experiments are being done mostly on layered ruthenates like , layered antiferromagnet like , layered manganites and some layered cuprates [1, 2, 3, 10, 11, 21]. The layered antiferromagnet, being spinPeierls compound, is proposed to exhibit a BerezinskiiKosterlitzThouless transition in the vicinity of spinPeierls transition temperature. The two dimensional spin half XY model was investigated and the validity of BerezinskiiKosterlitzThouless transition was confirmed [22]. The existence of Berezinskii Kosterlitz Thouless transition was proposed long ago in (S= layered ferromagnet) [6]. According to our knowledge layered ferromagnets with spin are the least studied systems, both from theoretical and experimental point of view, till date.
An extensive experimental study of spindynamics in a layered ferromagnet has been carried out by Hirakawa et. al. [6] using neutron scattering probe on . Their results exhibit a central peak (at ) in the plot of “neutron count vs. frequency” at a fixed value of the wavevector . Subsequent developments of approximate analytical theories and Monte Carlo Molecular Dynamics(MCMD) analysis have suggested that the existence of central peaks is partly due to scattering of neutrons from moving vortices and antivortices [14, 15].
Here we aim to examine how far the picture of ideal gas of vortices and antivortices could be extended to the quantum spin models. For this purpose, we choose as the reference system. This is a spin quasitwodimensional ferromagnetic material, on which extensive neutron scattering studies have been done.
The plan of the paper is as follows. In Sec. 2, we briefly describe the classical theory of mobile vortices and antivortices. In the same Section we explain our mathematical formulations in detail. In Sec. 3 we discuss our calculations and results. In Sec. 4 we present the conclusions and the future plan.
2 Mathematical Formulation
The dynamics of mobile vortices in a ferromagnetic system has already been treated both analytically and numerically by Huber [13] and Mertens et. al. [14]. In this work we apply their classical formalism to study the phase transition in . We calculate the spinspin correlations taking the experimental situations into account. For our purpose we present a brief description of the analytical treatment developed in Ref. [14]. The starting Hamiltonian is,
(1) 
where label the nearest neighbour sites on a two dimensional square lattice, J is the coupling constant and the classical spin vector is . This is an anisotropic Heisenberg Hamiltonian which, for , represents a ferromagnetic system. The quantity is the anisotropy parameter whose XY and isotropic Heisenberg limit correspond to and respectively. The general time dependent spin configuration in spherical polar coordinate system is given by,
(2) 
with . Following the formulation of Hikami and Tsuneto, the solutions are given by [20], and
(3)  
for single vortex centred at , where (3) describes the asymptotic behaviour of . Here vortex core radius is given by [15], . This type of spin configuration defines a ‘meronic’ type of the spin vortex.
The definition of the spinspin correlation function is given by,
(4)  
where represents the thermal average. In the case of classical gas of ideal vortices the thermal average has to be done by taking Maxwellian velocity distribution function. Here and are inplane correlations and is the outofplane correlation.
The effective analytical expression for the inplane correlation can be taken as [14],
(5) 
with , where is the root mean square velocity. Here is the vortex vortex correlation length. The rootmean squared velocity of the vortices was first calculated by Huber as [20];
(6) 
where is the density of free vortices at . The Fourier transform of in (5) gives rise to the inplane dynamical structure function given by,
(7) 
This is a squared Lorentzian, peaked at , with q dependent width,
(8) 
Exactly same results holds for also.
From the definition of , it can be shown that the outofplane correlation is given by [14],
(9) 
where is the Maxwell velocity distribution for a single vortex. Performing first the spatial Fourier transform and then the temporal Fourier transform, it can be shown that, the outofplane dynamical structure function has the form,
(10) 
Here is the velocity independent vortex form factor and it has the form . The form of as in (10) exhibits a central peak at . The width of the central peak is i.e., linear in .
In a typical inelastic neutron scattering experiment the count rate is related to the dynamical structure function as [23],
(11) 
At finite temperature there always exist creation and annihilation of excitations. A detailed balance condition is always needed to relate the intensities of up scattering () and down scattering (). True quantum mechanical , denoted by is recovered by the relation ,
(12) 
Where the factor is called the Windsor factor [17]. This incorporates the detailed balance condition, as required by the thermal equillibrium. Another important factor, which has to be taken into account, is the instrumental resolution function or . This essentially incorporates the different independent instrumental properties that affect the incident and scattered beam of neutrons [23]. In order to compare theory with experiment one has to convolute the theoretical expression, obtained from a model under consideration, by the resolution function [24]. Thus we consider the convoluted dynamical structure function given by,
(13)  
The resolution function has to be chosen so as to give minimum ripples at the end points of the resolution width. For this purpose a suitable window function such as Tukey window function may be chosen,
(14)  
The parameter , occurring in the Window function, can be set from the resolution half width obtained from experimental data.
2.1 InPlane Dynamical Structure Function
In our formulation for the inplane dynamical structure function we take into account the Tukey window function, as mentioned above (see equations (14)). Using (5), (13) and (14) we compute the Fourier transform of inplane spinspin correlation,
(15)  
Now, , where is Bessel function of order zero. The spatial integration is performed from zero to a certain radius . A final expression for the convoluted inplane dynamical structure function takes the form,
(16)  
Since, and are both even function in , only contributes to the temporal part of the integration. From symmetry Y component of the inplane dynamical structure function is same as X component of the inplane dynamical structure function . Let us note that in the above analysis the formulation holds only for . For the vortex vortex correlation length is not defined and hence the formalism can’t be extrapolated below .
2.2 OutofPlane Dynamical Structure Function
The outofplane dynamical structure function is given by,
(17) 
where is the Fourier transform of R(t). The reason for taking (17) as the expression for convoluted outofplane dynamical structure function is that, unlike (5), an analytical expression for can’t be evaluated from (9). So one has to start from (10).
The integral in (17) has been computed numerically. The defined here corresponds only to the mobile vortices; whereas the experimental data contain the contributions from bound vortices and fragile ‘spin wave like’ modes. These fragile ‘spin wave like’ modes are the largely decaying spin wave modes above the ferromagnetic paramagnetic transition temperature (Curie temperature). In order to compare with the experimental observations, one has to extract the mobile vortex contribution from the experimental data. This can be done by subtracting the fragile mode contribution and the frozen vortex contribution from the experimental data. The fragile mode contribution has been subtracted by taking the fragile mode contribution above transition temperature to be same as the spin wave contribution just below transition temperature. This is valid as long as we are considering the temperature which are not far below or above from the transition temperature. To find the approximate analytical expression for due to bound vortex contribution, the limiting value of is taken as in (10). Then from (10) it easy to find an expression for namely,
(18) 
Where is the bound vortex density. Since, the system has no net topological charge we can assume that there are equal number of vortices and antivortices present in the system and we can take assuming square lattice structure. This is correct as long as the temperature is just below where all the vortices are frozen but once the temperature crosses some of the bound vortices become mobile and the bound vortex density can be approximated as,
(19) 
where, is in the units of inverse of plaquette size(). Since, [14], given by (19) is temperature dependent. Here is of the order of lattice parameter. Using (18) and (19) the bound vortex contribution has to be subtracted carefully from the experimentally observed count.
We would like to point out that we could not apply this procedure to extract out bound vortex contributions in the case of inplane dynamical structure function (see Sec: 3).
2.3 Total Dynamical Structure Function (SpinSpin Correlation)
The general expression for the total dynamical structure function is,
(20) 
where, the total spinspin correlation is is defined by (4). So, the total dynamical structure function is . Since, X and Y components of the spins are symmetric we have, and the total dynamical structure function takes the form,
(21) 
Here, we would consider (21) only for mobile vortices.
It is an important fact that the formalism explained above incorporates the Windsor factor and the presence of in the quantum expression of magnetic moment corresponding to the spins constituting the vortex [13]. Therefore the formalism looks like a semiclassical one. Henceforth we will call our combined theoretical approach ‘semiclassical like’.
3 Calculations and Results
We apply the formalism of Sec: 2 on a real material for which neutron scattering experiments have been performed [6]. It is a quasitwodimensional spin ferromagnet, where the interaction is mainly Heisenberg type with only 1 XY like anisotropy . The transition is close to KT type with slight modification due to Heisenberg type interaction. The magnetic lattice structure for is approximately a body centred tetragonal lattice i.e. a lattice, composed of stacking of 2D square lattices [6]. The physical parameters are given in the Table 1, which have been used throughout the calculation.
parameter  magnitude 

exchange coupling (J)  11.93 K 
lattice parameter ()  4.123 Å 
‘b’  1.5 
‘’  5.5 K 
We start with the investigation of the inplane correlation (inplane dynamical structure function) . The radius in (16) is for a lattice, as used in the MCMD analysis by Mertens et. al. [14], where is the lattice parameter. We set the value of according to the experimental resolution width (0.01 meV) [6]. We compute numerically, for two different temperatures, 6.25 K and 6.75 K, for q(planar)= 0.04 reciprocal lattice units(in the units of ), experimentally being the momentum transfer.
There are two threshold values of [6], namely and , where for the system behaves like 2D Heisenberg system and for the system behaves as 3D XY system. For the system behaves as 2D XY system.We have varied the energy transfer from 0.3 meV to +0.3 meV, which includes the range 0.2 meV to +0.2 meV as taken in experiment [6]. The convoluted inplane dynamical structure function is plotted in Fig. 1 and Fig. 2, where is the natural time unit for the system/material (in our case ). These figures indicate that after convoluting with the Tukey window function, the inplane dynamical structure function no longer remains squared Lorentzian, though in both the cases central peaks persist. The width of the curve is much larger than that of the squared Lorentzian.
We notice that the convoluted inplane dynamical structure function function has become negative just above 0.1 meV. The occurence of negative values of the dynamical structure function has been dealt with in detail in sec.4 and in Appendix.
Again comparing Fig. 1 and Fig. 2 we find that the width of the squared Lorentzian increases with the increase of temperature whereas that of the does not undergo any change. Later, we will present a comparison of the convoluted total dynamical structure function with the experimental one(See Fig. 5 and Fig. 6).
We now evaluate the outofplane dynamical structure function for two different temperatures, 6.25 K and 6.75 K, for q(planar ) = 0.04 r.l.u, using (17). The expression for is,
(22) 
We use the same value of as used for . Here also the reasons for the choice of temperatures and q(planar) are same as that for the inplane dynamical structure function. In Fig. 3 and Fig. 4 we have plotted the outofplane correlation, . We have varied the from to in (17), where is estimated from the resolution width as before.
In this case the bound vortex contribution has been subtracted carefully, using (18) and (19), from the observed count at 6.25 K to obtain the effective mobile vortex contribution. The methodology for extracting the mobile vortex contributions from the experimental data has been explained in Sec 2.2. As long as the counts at 6.75 K are concerned, the fragile ‘spin wave like’ modes are highly decaying so that it can’t be assumed to be the same as the true spin wave modes observed at 5 K. So only bound vortex contribution has been subtracted at 6.75 K. The normalization factors, required for the quantitative comparison between the theoretical and the experimental results, have been estimated from the neutron count extracted from the experiment on [6].
We find that the outofplane dynamical structure function is also negative within the resolution width(see Figs. 3 and 4)! Moreover experimental peak is out side the resolution width, while the peak corresponding to the is at .
The above calculations lead us to the theoretical estimate for the convoluted total dynamical structure function given by (21). In Figs. 5 and 6, has been compared with the filtered experimental data obtained by subtracting the bound vortex contributions and fragile ‘spin wave like’ contributions (see Sec. 2.1 and 2.2). In these plots the intensities of the experimental peak and that of the central peak of the have been matched.
It is clear from Fig. 5 that at 6.25 K the experimental peak occurs approximately at 0.08 meV, which is way outside the resolution width. At 6.75 K [see Fig. 6] the peak of the experimental graph is not far from the central peak. It is reasonable to say that as the temperature is increased, we are getting better agreement of the with the experimental observations. This agreement isregarding the position of the central peak. Apart from the central peak there are two other peaks at finite frequency at both the temperatures. These are nothing but the reminiscent of the outofplane dynamical structure function contribution as seen from Figs. 3, 4, 5 and 6. This signifies the fact that the inplane correlation is largely dominating over the outofplane correlation.
The total spinspin correlation is still negative just above 0.1 meV.
Though it is true that the dynamical structure function can’t be negative, here in our case the negativity occurs as a result of the convolution of analytical expression of . Even for a conventional long range ordered system, the dynamical structure function corresponding to a classical pure spin wave comes out to be negative beyond a certain range of frequency when convoluted with any spectral window function. Furthermore the above peculerity persists even when quantum effects are incorporated through a detailed balance factor[see Appendix].
The inclusion of quantum mechanical detailed balance factor in the semiclassical like treatment for dynamics of mobile vortices and antivortices, is not even causing any appreciable asymmetry, as seen in the theoretical plots in our case of spin. The theoretical plots are largely symmetric around . A very small asymmetry in the theoretical plots are being seen for higher values of while the experimental data are showing clearly the asymmetry.
It may be noted that in our analysis the bound vortex contributions have been approximately estimated only for outofplane dynamical structure function . This is because, in this case, we are able to truncate the expression, as given in (10), to the regime , by making . In (10), there exists no explicit dependence of on the correlation length . In case of inplane correlation, as given in (5), we need to find for due to its explicit appearance in that expression. Since is not defined for we are not able to estimate the bound vortex contribution for inplane dynamical structure function.
In summary, we find that the width of the convoluted inplane dynamical structure function is much larger than that of the squared Lorentzian. Values of the inplane dynamical structure function comes out to be negative beyond a finite range of energy transfer. The convoluted outofplane dynamical structure function becomes negative as well; however this happens within the resolution width about the central peak (peak at ). The total convoluted dynamical structure function also becomes negative in the regime where the inplane dynamical structure function had become negative. No appreciable asymmetry is created even after including the Windsor factor. We find that for both the temperatures the convoluted total dynamical structure function is symmetric around ; whereas the experimental observation is not. The theoretical model of semiclassical treatment of ideal gas of unbound vortices tends to agree with the experimental observations better at higher temperatures (for spin system), when we consider the experimental results at T= 6.25 K and T= 6.75 K (Figs. 5 and 6). It is worthwhile to point out that same results hold for unbound antivortices also.
4 Conclusions & Discussions
The laws of quantum mechanics which govern all real systems, ensure the dynamical structure functions to be always positive definite [23]. We find in our analysis however, that the semiclassical treatment based on ideal gas of vortices (antivortices) for a low spin system leads to the occurence of negative values of dynamical structure function, over a large range of energy transfer, when convoluted with any standerd resolution function.
Based on the analysis carried out in the Appendix we can infer that for the dynamics of mobile vortices and antivortices, the negative values of are occurring due to the following factors;
i) the choice of the resolution function, whicn in our case is the Tukey function,
ii) the choice of the value of resolution width , which in this case is made fixed by experimentally imposed resolution width,
iii) use of a semiclassical like treatment to extend the classical theory of dynamics of mobile vortices (antivortices) to a quasi two dimensional spin ferromagnet which is quantum mechanical.
To avoid the negativity in the we could have chosen a different resolution function. Indeed, it has been shown that most of the resolution functions are more or less oscillatory in the fourier space [24, 29] . An extra smoothening factor can be used to dampen the oscillation of the resolution function. This extra factor is eventually related to the resolution width and it makes the resolution function smoother if the resolution width is decreased [24]. However, in our case the resolution width is fixed from the experiment and consequently the oscillation of the resolution function can’t be avoided by merely changing the resolution function.
It has been found that there is a range of over which remains positive. We can call it as the physically admissible range. This range is related to the magnitude of the spin occuring in the theoretical model under consideration and to the resolution width. On the basis of the analysis (presented in the Appendix) it is expected that even in the case of dynamics of mobile vortices and antivortices the physically admissible range would be larger for higher spin value and smaller for lower spin values.
Another way to avoid the negativity is to assume, outside the physically admissible range [24]. If this physically admissible range is within the range of experimental interest then the assumption is not applicable. In our case of spinhalf ferromagnet the physically admissible range is well within the range of over which the neutron scattering data has been taken in the experiment (as seen from Figs. 5 and 6).
Hence, the negative values of can only be due to the use of the semiclassical like treatment to extend the classical theory of dynamics of mobile vortices and antivortices to a quasi two dimensional spin ferromagnet ().
Moreover, the convoluted outofplane dynamical structure function computed from our semiclassical like treatment, becomes negative within the experimentally imposed resolution width itself. Thus the central peak occurring in this case may not posses a well defined width.
The agreement between the behaviour of dynamical structure functions obtained from our theoretical calculations and that from the experiment, in terms of the peak position and the overall shape, is found to be fairly good at temperatures much larger than .
Although the vortices (and antivortices) are extended objects the Maxwell Boltzmann distribution can still be used for the motion of the centre of mass of these objects.
Our investigation presented in this paper brings out the fact that a complete quantum mechanical treatment is essential for describing the detailed features of the dynamics of unbound spin vortices and antivortices corresponding to low spin magnetism systems. As a first step towards this, a theoretical framework for describing static quantum spin vortices and antivortices and their topological properties has been developed [25, 26, 27, 28]. An extension of this formalism to the case of mobile spin vortices and antivortices is crucial for the quantum mechanical calculation of dynamical structure function. This would go a long way towards an explanation for the experimental results observed for the genuine quantum spin systems like .
5 Acknowledgments
One of the authors (SS) acknowledges the financial support through Junior Research Fellowship (09/ 575 (0089) / 2010 EMR–1) provided by Council of Scientific & Industrial Research (CSIR) and also acknowledges Mr. Bandan Chakrabortty for valuable discussions regarding the computer programming.
Appendix: corresponding to classical spinwave
For classical spin wave at very low temperature corresponding to a classical Heisenberg ferromagnet, the dynamical structure function has the form:
(23)  
Where, for a cubic lattice. Here is the number of nearest neighbours and .
Convoluted dynamical structure function, has been defined in (13), where is the fourier transform of a siutably chosen spectral function. So far Tukey function has been used and the same is used in this case also. The fourier transform of Tukey function is given by (22). Using (13), (22) and (23) we find the convoluted dynamical structure function corresponding to spin wave, to be,
(24) 
Analyzing (24) we find that there exists a range
(25) 
within which remains positive and outside, it becomes negative. This remains true even if we choose any other resolution function. This is basically due to the fact that when fourier transform is performed on spectral functions (defined in time domain), the resulting functions in space mostly turn out to be oscillatory [29].
One way to avoid these negative values of would be to assume outside [24]. This prescription however, can’t be taken into consideration when comparing theoretical predictions with experimental results, if the energy range of experimental interest contains the above mentioned range of . Another way to avoid the negative values of is to decrease the value of [see (25)], where is related to the experimental resolution width () by the relation,
(26) 
being in energy units. Decrease in means increase in the value of which signifies a poor experimental resolution width. In comparing the theoretical predictions with the experimental results, is determined by the relation (26), where is fixed from the experiment.
Another interesting feature of (25) is the fact that the region of where remains positive, depends on the spin value S via the relation,
(27) 
This shows that when is fixed from the experimental resolution width, S remains as the free parameter to determine the range of as mentioned in (25). Increase in S will increase the range of where remains positive. Hence for a system with high spin value S, the range of positive may actually be at par with the energy range of experimental interest, unlike the case of spin .
It can be shown from (24) that the range of , given by (25), will not change even if a detailed balance factor (in this case Windsor factor) is introduced to incorporate the quantum effects.
References
 [1] A. A. Katanin, V. Yu Irkhin, Physcis– Uspekhi, 50(6), 613 635(2007) and references therein.
 [2] J. Prokop et. al., Phys. Rev. Lett 102 177206(2009), M. Cwik et. al., Phys. Rev. Lett 102 057201(2009); F. Demmel and T. Chatterji, Phys. Rev. B 76 212402(2007); T. Chatterji et. al., Phys. Rev. B 76 144406(2007).

[3]
F. Weber et. al., Phys. Rev. Lett 107 207202(2011); S. M. Yusuf et. al., Phys. Rev. B 84 064407(2011); H. Ulbrich et. al., Phys. Rev. B 84 094453(2011); H. J. Lewtas et. al., Phys. Rev. B 82 184420(2010); L. Braicovich et. al., Phys. Rev. B 81 174533(2010);
Ranjan Chaudhury, Spin dynamics of quantum spin models in layered systems– signature of Kosterlitz– Thouless scenario? – presented at ULT 2011 (August 19– 22, 2011), held at KAIST, Daejeon, South Korea;
Yong– il Shin, Creation of 2D Skyrmions in spin– 1 polar Bose– Einstien Condensates – presented at ULT 2011 (August 19– 22, 2011), held at KAIST, Daejeon, South Korea.  [4] J. K. Kjemsand and M. Steiner, Phys. Rev. Lett., 41, 1137(1978); I. U. Heilmann, et.al., Phys. Rev. B., 24, 3939(1981); G. Reiter, Phys. Rev. Lett., 46, 202(1981).
 [5] H. Bethe, Z. Physik 71, 205(1931); Y. Endoh, et.al., Phys. Rev. Lett., 32, 170(1974); F. D. M. Haldane, Phys. Rev. Lett., 50, 1153(1983).
 [6] K. Hirakawa, H. Yoshizawa, K. Ubukoshi, J. Phys. Soc. Jpn., 51, 2151(1982); K. Hirakawa, et.al., J. Phys. Soc. Jpn. Part2, 52, 4220(1983).
 [7] M. T. Hutchings, J. AlsNielsen, P. A. Lindgard and P. J. Walker, J. Phys. C: Solid State Phys., 14, 53275345(1981).
 [8] M. T. Hutchings, P. Day, E. Janke and R. Pynn, J. Magn. Magn. Mater, 5457, 673(1986).
 [9] D. G. Wiesler, et.al., Physica 136B, 2224(1986).
 [10] P. Steffens et.al., Phys. Rev. B., 83, 054429(2011).
 [11] L. Capogna et. al.; Phys. Rev. B., 67, 012504(2003).
 [12] K. Yamada et.al., Phys. Rev. B., 40, 4557(1989).
 [13] D. L. Huber, Phys. Lett. 68A, 125(1978); Phys. Rev. B. 26, 3758(1982).
 [14] F. G Mertens, A. R. Bishop, G. M. Wysin and C. Kawabata, Phys. Rev. Lett, 59, 117(1987); Phys. Rev. B. 39, 591(1989).
 [15] M. E. Gouva, G. M. Wysin, A. R. Bishop, F. G. Mertens, Phys. Rev. B, 39, 11840(1989).
 [16] A. R. Volkel, G. M. Wysin, A. R. Bishop and F. G. Mertens, Phys. Rev. B. 44, 10066(1991).
 [17] M. T. Evans and C. G. Windsor, J. Phys. C., 6, 495 (1973); C. G. Windsor abd J. LockeWheaton, J. Phys. C., 9, 2749 (1973); R. Chaudhury and B. S. Shastry, Phys. Rev. B. 37, 5216(1988).
 [18] J. M. Kosterlitz and D. J. Thouless: J. Phys. C 6 1181(1973).
 [19] V. L. Berezinskii, Sov. Phys. JETP 32, 493(1970); Sov. Phys. JETP 34, 610(1972).
 [20] S. Hikami and T. Tsuneto, Prog. Theor. Phys. 63, 387(1980).
 [21] J. E. Lorenzo et. al., Europhys. Lett., 45(1), 4551(1999).
 [22] K. Harada and N. Kawashima, J. Phys. Soc. Jpn., 67, 2768(1998).
 [23] S. W. Lovesey, Theory of Neutron Scattering from Condensed Matter(Oxford– Clarendon Press) 1986; F. Hippert, E. Geissler, Jean. L. Hodeau, E. LBerna and JeanR Regnard, Neutron and XRay Spectroscopy, Grenoble Sciences, Springer 2006.
 [24] M. Takahashi, J. Phys. Soc. Jpn., 52, 3592(1983).
 [25] E. Fradkin, M. Stone, Phys. Rev. B., 38, 7215, (1988); E. Fradkin, Field Theoies of Condensed Matter Systems, (AddisonWesley, CA, 1991).
 [26] F. D. M. Halden, Phys. Rev. Lett, 61, 1029(1988).
 [27] E. Jr. Loh, D. J. Scalapino, P. M. Grant, Phys. Rev. B., 31,4712(1985); D. D. betts, F. C. Salevsky, J. Rogiers, J. Phys. A, 14, 531(1981).
 [28] R. Chaudhury and S. K. Paul, Phys. Rev. B., 60, 6234(1999); Mod. Phys. Lett. B, 16, 251259(2002); Eur. Phys. J. B, 76, 391398(2010).
 [29] F. J. Harris, Proceedings of The IEEE, 66, 51(1978).