{'fraction_non_alphanumeric': 0.0488775706734576, 'fraction_numerical': 0.01697018981788611, 'mean_word_length': 4.063935847420893, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Critical dynamics in the 2d classical XY-model: a spin dynamics study

Using spin-dynamics techniques we have performed large-scale computer simulations of the dynamic behavior of the classical three component XY-model (i.e. the anisotropic limit of an easy-plane Heisenberg ferromagnet), on square lattices of size up to 192^2, for several temperatures below, at, and above T_KT. The temporal evolution of spin configurations was determined numerically from coupled equations of motion for individual spins by a fourth order predictor-corrector method, with initial spin configurations generated by a hybrid Monte Carlo algorithm. The neutron scattering function S(q,omega) was calculated from the resultant space-time displaced spin-spin correlation function. Pronounced spin-wave peaks were found both in the in-plane and the out-of-plane scattering function over a wide range of temperatures. The in-plane scattering function S^xx also has a large number of clear but weak additional peaks, which we interpret to come from two-spin-wave scattering. In addition, we observed a small central peak in S^xx, even at temperatures well below the phase transition. We used dynamic finite size scaling theory to extract the dynamic critical exponent z. We find z=1.00(4) for all T<= T_KT, in excellent agreement with theoretical predictions, although the shape of S(q,omega) is not well described by current theory.


Introduction
The two-dimensional XY -model is one of the 'special' models of magnetism. It undergoes an unusual phase transition to a state with bound, topological excitations (vortex pairs) but no long range order [1]. The (three component) XY-model may be viewed as a special case of the anisotropic Heisenberg model in which the coupling between the zcomponents of spins vanishes. It has static properties which are similar to those of the "plane rotator" model, in which the spins have only two components. The static properties which of both models have been determined by numerical simulation [2,3,4] and found to be consistent with the predictions of the Kosterlitz-Thouless theory. For example, the susceptibility shows an essential singularity instead of a power law divergence, and computer simulations show that vortex pairs unbind at T KT . The model is critical, i.e. it has infinite correlation length ξ, at all temperatures T ≤ T KT . The spin-spin correlation function decays algebraically with distance for all T ≤ T KT , but with a power η which varies with temperature.
The dynamic behavior of the model should be governed by the dynamic critical exponent z, which describes the divergent behavior of the relevant time scale [5], i.e. τ ∝ ξ z for all T ≤ T KT . Recently, finite size scaling for critical dynamics in the neutron scattering function has been developed [6], and successfully applied to the study of a 3-dim. Heisenberg ferromagnet. The XY -model has true dynamics which can be determined by integrating the equations of motion for each spin; its critical dynamics has been studied theoretically [7,8,9,10] with different predictions for the nature of the neutron scattering function. In contrast, the "plane rotator" model does not possess equations of motion; it has only stochastic, i.e. purely relaxational time dependence, which has also been examined by Monte Carlo simulation [4]. Note, however, that different dynamic exponents are expected for stochastic and for true dynamics [5].
In this paper we present the first large scale simulation study of the true dynamic behavior of the XY -model. Great care was taken to ensure that statistical as well as systematic errors were both well understood and quite small. An earlier, much less complete study [11] indicated a rich structure in the neutron scattering function which was not adequately described by theory. Our model is defined by the Hamiltonian where S i is a three-component classical vector of length unity and the sum is over all nearest neighbor pairs. The equation of motion of each spin is [11] d dt where andê x andê y are unit vectors in the x-and y-directions respectively. Eq. (2) represents a set of coupled equations and can be integrated numerically. (The plane rotator model has a hamiltonian of the same form as eq. (1), but since the vectors have only two components, an equation of motion cannot be defined in the same way as for the XY -model.) In section 2 we define the neutron scattering function, provide dynamic finite size scaling equations, and summarize analytical results. Section 3 describes the details of our simulations. We present our data and subsequent analysis in section 4, and draw conclusions in section 5.
The results of our dynamic study prompted us to reexamine the static properties of the model. In order to obtain more reliable values for the critical temperature and the static critical exponent η, we carried out static Monte Carlo simulations; our results are presented in the appendix.

Neutron Scattering Function S(q, ω)
The neutron scattering function S(q, ω) (also called the dynamic structure factor) is an experimental observable and is fundamental to the study of spin dynamics. It is defined [6] for momentum transfer q and frequency transfer ω as the space-time fourier transform of the space-displaced, time-displaced spin-spin correlation function where k = x, y, or z is the spin component, displacement r is in units of lattice spacings, and the angle brackets ... denote the thermal ensemble average. Note that in the 2-dim. XY -model, S k r (t) ≡ 0 for all components k = x, y, z. The equations of motion (2) are time reversal invariant, therefore C kk (r − r ′ , t − t ′ ) is symmetric in t and t ′ , and S(q, ω) is real-valued.
The neutron scattering function generally depends on the correlation length ξ and may be written in the form [12] where S kk ξ (q) is the total intensity given by and f kk is a normalized shape function, ∞ −∞ f kk (x, q, ξ)dx = 1. The characteristic frequency ω kk m (q, ξ) is a median frequency determined by the constraint In dynamic scaling theory it is assumed that the median frequency ω kk m (q, ξ) is a homogeneous function of q and ξ, i.e., where z kk is the dynamic critical exponent and Ω kk is another shape function, and that the function f kk depends only on the product of qξ but not on q and ξ separately. Therefore S kk ξ (q, ω) as given in eq. (6) simplifies to Note that the z-component of the magnetization is conserved during the time evolution. Thus the neutron scattering function S kk ξ (q, ω) can be regrouped in terms of symmetry into the out-of-plane component S zz and the in-plane component with different physical behavior. (As reported in sections 4.3 and 4.4, we find that for the two different components the exponent z is the same, but the scaling functions Ω and f differ.)

Dynamic Finite Size Scaling
At the critical temperature T KT and below, the XY -model is expected to be critical, with infinite correlation length ξ. In this region the dynamic critical exponent z can be extracted by using the dynamic finite size scaling theory developed by Chen and Landau [6]. These authors also introduced a frequency resolution function in order to smoothen the effects of finite length of time integration, which was not necessary for most of the analyses in our study because of much longer integration times (see section 3). Their dynamic finite size scaling relations [6] can then be simplified to (replacing a factor ω in front of S kk L (q, ω)), and analogous to eq. (9) and eq. (10). We see that two different ways emerge to test dynamic scaling and to estimate the dynamic critical exponent z : Firstly, from eq. (13), z is given by the slope of a graph of log ω m versus log L at fixed value of qL. Secondly, eq. (12) implies that for correctly chosen z and at fixed value of qL, graphs of (S kk L (q, ω) / L z S kk L (q)) versus ωL z should all fall onto the same curve for different lattice sizes. Both procedures will only be valid for sufficiently large lattice size.

Analytical Results
The dynamics of two-dimensional ferromagnets with easy-plane asymmetry, specifically the XY -model, were first analyzed by Villain [7] and by Moussa and Villain [13]. The in-plane scattering function was found to have a delta-function spin-wave peak at low temperature, and a spin-wave peak of the form close to T KT . Here η(T ) is the critical exponent describing the decay of the static spin-spin correlation function (and we now expect η = 1 4 at T KT ). Nelson and Fisher [8] treated the classical XY model in a fixed length hydrodynamic description for T ≤ T KT , without vortex contributions. They obtained the transverse spin-spin correlation function, for r, ct ≫ 1 and ct = r, where c is the spin-wave velocity. The fourier transform of eq. (15) has the form where the scaling function Ψ behaves like around the spin-wave peak, and for large values of ω/cq. Nelson and Fisher also predicted that the dynamic critical exponent z, eq. (9), which is expected [5] to be z = d/2 for d > 2, is Note that the value z = 1 and a linear dispersion relation are also implicit in the argument ω/cq of the scaling function Ψ in eq. (16). Finally, both Villain and Nelson and Fisher predicted a very narrow (delta-function) spin-wave peak in the out-of-plane function S zz , at ω = cq.
More recently, Menezes, Pires, and Gouvêa [9] have performed a low temperature calculation which includes the contribution of out-of-plane fluctuations to the in-plane correlation functions. They worked in the harmonic spin-wave approximation which is justifiable for large spin s and therefore also for our classical spins, and they used a projection operator technique. They found a spin-wave peak similar to that of Nelson and Fisher, whereω = ω/(cq). In addition to the spin-wave peak, they found a logarithmically diverging central peak, i.e. a signal at very small ω, which diverges like Of course, a central peak at small temperature can also be caused by other mechanisms, for example vortex pairs diffusing like a dilute pair of solitons [10].
The dynamic behavior of the XY -model is different above T KT . For a phase transition of Kosterlitz-Thouless type, the spin stiffness should drop discontinuously to zero at T KT , i.e. the spin-wave peak is predicted to disappear [8,14]. Above T KT , vortex-antivortex pairs unbind, and their diffusion leads to a strong central peak in S(q, ω).
Mertens et al. [15,16] calculated S(q, ω) above T KT , assuming an ideal dilute gas of unbound vortices moving in the presence of renormalized spin-waves, and screened by the remaining vortex-antivortex pairs. They found a Lorentzian central peak for S xx , and a Gaussian central peak for S zz , where γ = 1 2 √ πū/ξ,ū is the rms vortex velocity, ξ the correlation length, and n v ∼ (2ξ) −2 the free vortex density; and they compared their results to numerical simulations (see below).

Previous numerical work
Gerling and Landau [11] carried out spin dynamics simulations on the XY model with L ≤ 204 and found both spin-wave peaks and a central peak. The resolution was too limited, however, to allow quantitative comparison with theory or to extract an estimate for the dynamic exponent. Mertens et al. [15,16] performed spin dynamics simulations, with fixed system size L = 100 and very low statistics (3 starting configurations) at T = 0.5 and T = 1.1. Below T c they observed only spin-wave peaks in both S xx and S zz ; above T c they saw a strong central peak and a weak spin-wave peak in S xx , and vice versa in S zz . The width and intensity of the central peaks were compatible with eq. (22) and eq. (23).
Other earlier numerical work on dynamical behavior has been exploratory [17,18].

Experiments
The closest physical realizations of the XY-model are materials with very large anisotropy, more generally described by strongly anisotropic Heisenberg models. Several experiments have studied the dynamics of such materials [19,20], like Rb 2 CrCl 4 , K 2 CuF 4 , and CoCl 2 .
In a recent study on stage-2 CoCl 2 intercalated graphite, Wiesler et al. [19] found four temperature regimes with different behavior. There are indications of a Kosterlitz-Thouless transition at a temperature "T u ", though some properties disagree with KT predictions. Between temperatures "T l " and "T u ", they observed spin-wave peaks. It is not clear whether ta central peak is present there. (In this temperature region the long range part of the scattering function shows true 2-dimensional character, whereas for T < T l 3-dimensional correlations develop.) Above "T u ", the in-plane scattering function showed the expected central peak, and the out-of-plane function exhibited damped spinwaves.
In experimentally available materials both defects and the effects of residual threedimensional couplings limit the effective size of the two-dimensional KT-like system to L exp = O(100) lattice spacings [21,19]. Remarkably, this size is similar to the lattice sizes of the present numerical study. For further discussion and an extensive listing of relevant literature, see the recent overview contained in [19].

Simulations
We have studied the two-dimensional classical XY -model with Hamiltonian given in eq.
(1) on L×L lattices with periodic boundary conditions for 16 ≤ L ≤ 192, at temperatures T = 0.4, 0.6, 0.7, 0.725, and 0.8 in units of J/k B . Most of these temperatures are in the critical region T < ∼ T KT = 0.700(5) (see appendix). Equilibrium configurations were created at each temperature using a Monte Carlo method which combined cluster updates of the x and y spin components (using the Wolff embedding method [22,23]) with vectorized Metropolis and overrelaxation [24] spin reorientations. After each single-cluster update, two Metropolis and eight overrelaxation sweeps were performed. Use of the cluster algorithm was important, since critical slowing down was severe for most of our simulations; the inclusion of cluster flipping reduced Monte Carlo autocorrelation times at L = 192 and T = 0.6 from more than 300 to about 3 hybrid sweeps, while requiring only a factor of two more CPU time per sweep. We performed 200 hybrid sweeps between equilibrium configurations, and discarded the first 5000 hybrid sweeps for equilibration.
Between 500 and 1200 equilibrium configurations were generated for each lattice size and temperature. We found this many configurations to be necessary in order to sufficiently reduce statistical errors in the resulting neutron scattering function. The error bars in our figures represent the statistical errors for averages over the equilibrium configurations, drawn from the canonical ensemble.
Starting with each equilibrium configuration, the time dependence of the spins was determined from the coupled set of equations of motion, eq. (2), and was integrated numerically using a vectorized fourth order predictor-corrector method [6], with a time step size of δt = 0.01J −1 . The maximum integration time was generally t max = 400J −1 ; a few runs were also performed for lattice size 256 × 256 with t max = 800J −1 and produced the same physical results.
The time-displaced, space-displaced spin-spin correlation functions C(r − r ′ , t − t ′ ), eq. (5), were measured for each time integration, with and were then averaged. By fourier transformation in space and in time, eq. (4), we obtained the neutron scattering function S(q, ω). The time integration in eq. (4) was performed using Simpson's rule, with a time step of 0.1J, which has been shown [6] to be sufficiently small.
To reduce memory and computer time requirements, we restricted ourselves to momenta q = (q, 0) and (0, q), with q determined by the periodic boundary conditions, and data from these two spatially equivalent directions were averaged together to further enhance the statistical accuracy. We used fast fourier transforms to increase the efficiency of the program in calculating correlation functions. The frequency resolution ∆ω of our results is determined by the time integration cutoff t cutof f = 0.9 t max , see eq. (24), which will introduce oscillations of period 2π/t cutof f into S(q, ω). Since we observed very sharp spin-wave peaks (see section 4.1), we chose to integrate the equations of motion to very large times. We used t max = 200J −1 for L ≤ 96, and t max = 400J −1 for L ≥ 128. A theoretical delta-function in frequency will then become a widened peak with a width at half maximum of in the simulation data. To smoothen the oscillations, previous spin dynamics studies [11,6] have employed a frequency resolution function, replacing to compute S(q, ω). Because of the large values of t max in our study, we achieved a very small frequency resolution ∆ω, and the oscillations were not noticeable for most of our data. We therefore did not generally use a frequency resolution function and could significantly simplify the analysis of our data.
The numerical integration of the equations of motion can potentially become unstable at very large integration times. We checked that for our calculations, in which we integrate to much larger times than previous studies, we do not encounter this problem. We verified that the constants of motion (energy and magnetization in z-direction) do remain constant, with a relative variation of less than 3 × 10 −6 . We also verified that the neutron scattering function remains virtually unchanged when an additional integration of length t = 200J −1 is performed from each equilibrium configuration before starting to calculate time-displaced spin-spin correlation functions. All simulations were carried out using highly vectorized programs on the Cray C90 at the Pittsburgh Supercomputing Center.

Results
We now present our results for the dynamic structure factor S(q, ω), its dependence on temperature, frequency, momentum, and lattice size, and we analyze its dynamic scaling behavior. With few exceptions we have analyzed the data without the use of a frequency resolution function eq. (27). (The effects of such a function, and of integrating to shorter maximum times are described together with initial results in reference [25].) In order to investigate the XY -model in the critical phase, we chose several temperatures T ≤ T KT , including the best previous estimate of T KT ≈ 0.725J/k B [2], and one temperature well above the transition, T = 0.8J/k B . The results of our analysis prompted us to perform additional static Monte Carlo studies, which are described in the appendix. They provided an improved estimate of T KT = 0.700(5)J/k B . In order to elucidate the situation closer to the transition, we performed additional (but less extensive) spin dynamics simulations at T = 0.700J/k B . Figure 1 shows the temperature dependence of S(q, ω) as a function of ω, for lattice size L = 128 and fixed small momentum q = π/48 (i.e. n q = 2 in eq. (25)). Here, as in other results which we shall show, the error bars are determined from the statistical variation of results obtained from different initial spin states. The in-plane component S xx , fig. 1(a), exhibits a very strong and moderately sharp spin-wave peak at temperatures T < ∼ T KT . Even at the lowest temperature, however, the width of the peak is larger than the minimum value eq. (26) due to finite cutoff time. The position of the peak moves towards lower ω as the temperature increases, and the peak broadens slightly. Just above the transition, at T = 0.725, there is still both a strong spin-wave peak and a sizeable central peak in S xx . At higher temperature, the spin-wave peak disappears completely (for this low momentum) and only a large central peak remains. Note that from KT-theory [1] one would expect complete disappearance of a spin-wave peak at all T > T KT .

Spin-Wave Peak
There is sizeable additional structure in S xx away from the spin-wave peak at temper- atures up to T KT . We will discuss this structure in the following sub-section. The out-of-plane component S zz , shown in fig. 1(b), is two orders of magnitude weaker than the in-plane component. It exhibits a very sharp spin-wave peak for T ≤ T KT , whose width is limited by our ω-resolution. The finite time cutoff eq. (26) produces very noticeable oscillations in S(q, ω). (The magnitude of these oscillations is minute compared to the intensity of the spin-wave peaks in S xx .) The oscillations can be smoothened by convoluting S(q, ω) with a gaussian resolution function in frequency, as is shown in the inset. No central peaks are visible in S zz at T ≤ T KT . At T = 0.725, the peak in S zz is still present, with a larger width similar to that in S xx . In contrast to S xx , there is a clear, but weak, spin-wave peak in S zz even at T = 0.8 > T KT and small momentum. It is of similar intensity as the peak at lower temperatures. (See also section 2.5). Figure 2 shows the lattice size dependence of S(q, ω), at fixed momentum q = π/16. Below the transition, Fig. 2(a), the intensity of the spin-wave peak depends strongly on lattice size, whereas its position is constant. (The out-of-plane component S zz is dominated by finite time cutoff effects for T < T KT , and we do not show it here). Just above the transition, at T = 0.725, the spin-wave peak in S xx appears to gain intensity slightly as L increases, whereas neither its central peak nor the spin-wave peak in S zz show any finite size effects. At higher temperature, fig. 2(d) and 2(e), there is no visible lattice size dependence in either S xx or S zz . Notice the two different vertical scales for S xx and S zz . Data taken for L = 16 and L = 32 exhibit such strong finite size rounding that we have chosen not to show the data here.
In figure 3 we show the momentum dependence of S(q, ω). Fig. 3(a) and (b) display the behavior at T = 0.725, which is qualitatively similar to that at lower temperatures. The position of the spin-wave peak is the same for S xx and S zz and is proportional to momentum for small q. As q increases, the peak broadens, and becomes less intense, yet it remains well defined. The additional structure in S xx is strongly momentum-dependent, as will be discussed below.
For the zz-component, both the total intensity and the relative loss of intensity with increasing momentum are much smaller. Our ω-resolution dominates the width of the spin-wave peak in S zz only at the smallest q (which also appears in figure 1); it is not dominant at higher momenta or for S xx . We conclude that S zz has the expected deltafunction form only for very small momentum. (Higher order perturbation theory also predicts a finite width [26].) As shown in fig. 3(c), the intensity of the spin-wave peak decreases much more rapidly for S xx than for S zz with increasing q; the intensities cross each other well before the zone edge is reached. This behavior is similar at other temperatures. We also note that at all temperatures the total intensity S zz (q) (not shown) is constant with q, whereas S xx (q) decreases (as q η−2 ) and crosses S zz (q) at a slightly larger momentum q ≈ 1.5 − 2.
Well above the transition, at T = 0.8, fig. 3(d), S xx has no noticeable spin-wave peak at small momentum. The strong central peak rapidly loses intensity with increasing momentum. In marked contrast, the behavior of S zz ( fig. 3(e)) is very similar to that at lower temperatures, with clear but broadened spin-wave peaks. Notice that there is   now non-zero intensity at small ω in S zz . Remarkably, at very large momenta spin-waves appear in S xx even for T = 0.8 ( fig. 3(f), so that both a central peak and a spin-wave peak are present. Note that the vertical scale in fig. 3(f) is 100 times smaller than in fig.  3(d). There is no noticeable lattice size dependence here.  Figure 4 shows the position ω p of the spin-wave peak as a function of momentum. The expected linear portion of the dispersion curve extends to rather large momenta. With increasing temperature, the spin-wave velocity ∂ω p /∂q, which is proportional to the spin-wave stiffness, decreases slowly and approximately linearly, as shown in the inset, and theoretically expected for small T [14,27]. At T ≤ T KT , ω p is the same for S xx and for S zz , as expected by theory [7,8]. At T = 0.8 on the other hand, we can only plot the position of the residual peak in S zz , because S xx has dropped sharply to zero here, as expected for a KT transition.

Additional Structure in S xx (q, ω)
If we expand the vertical scale in plots of S xx (q, ω), we find that the in-plane component S xx (q, ω) shows rich structure in addition to the spin-wave peak. Note that the intensity of this structure is typically 10 −2 of the maximum. It is visible at all temperatures T < ∼ T KT . At the lowest temperature, the absolute intensity of this structure is low, but the relative intensity is quite high (see also fig. 11(a)). At higher temperature, the structure becomes rather smeared. No such structure can be found in S zz (q, ω).
The locations of additional maxima in S xx are essentially unchanged with L when n q = qL/(2π) is held fixed. Figures 5(b) and (c) show, on a logarithmic scale, that for odd values of n q there are strikingly regular pronounced peaks at ω = ω p /n q , 3ω p /n q , ..., and for even n q such peaks appear at ω = 0, 2ω p /n q , 4ω p /n q , .... At large n q , fig. 5(d), individual peaks cannot be distinguished; instead S xx is nearly constant below the spin-wave peak there. (In figs. 5(b) and (c), the data have been smoothened slightly with a resolution function eq. (27) with δω = 0.01, in order to reduce the wiggles and allow general features to be identified.) In addition to the regularly spaced pronounced peaks, there is further "fine structure" in S xx , clearly visible in fig. 5(a). In the course of our study, the additional structure became clearer as the statistical quality of the data improved; it is apparent that the structure is not statistical noise. Very close to the spin-wave peak, part of the additional structure may be due to the finite-time cutoff in our time-integration; but most of the observed structure must be due to different reasons.
One simple explanation for the observed rich structure, which is consistent with the data but for which we have no rigorous theory, is that of multi-spin-wave effects. Of these, two-spin-wave processes are likely to be the most important. Thus, at a given total momentum q we can have either a single spin-wave excitation of momentum q, or two spin-waves for which the sum or difference of momenta equals q. The result will then be both a single-spin-wave peak at a characteristic frequency ω p (q) as well as additional sum and difference peaks due to the two-spin-wave processes.
Of particular interest is the case when the two spin-waves have momenta q 1 and q 2 that are collinear, so that q = q 1 + q 2 is a scalar equation. Since the momenta are discrete on a finite lattice, q i = n q i 2π L , this implies n q = n q 1 + n q 2 . With a linear dispersion relation ω = cq, the difference of the two spin-wave frequencies is then ω = (2n q 1 − n q ) c 2π L , i.e. just the series of additional peaks that are visible in figures 5(b) and (c).
Using the measured dispersion relation ( fig. 4), we have calculated the frequencies of two-spin-wave excitations consisting of the most likely individual spin-waves, i.e. those with smallest individual momenta. For the case of n q = 3 and T = 0.6, these locations are marked in fig. 5(e). They are identified by the coordinates n q 1 of one of the spin-waves in reciprocal space; the sum of the two spin-wave momenta must equal q = (3 × 2π L , 0). The locations of the resultant excitations agree extremely well with the positions of the small peaks in S(q, ω), but we have no way of comparing intensities.
The presence of distinct small peaks at ω = 0, 2ω p /n q , 4ω p /n q , .... for even n q at T < T KT complicates the identification of a possible central peak. Interpolating the intensities for odd values of n q (which do not show peaks at ω = 0) to obtain estimates for even n q , we conclude that there is indeed extra intensity at ω = 0 which is not attributable to two-spin-wave processes.

Finite size scaling of characteristic frequency ω m
Equation (8) defines the characteristic frequency ω m of the whole spectrum of S(q, ω).
When there is only a single spin-wave peak, then ω m coincides with the spin-wave frequency ω p . This is the case at T = 0.4, where all frequencies coincide (within error bars), ω xx p = ω xx m = ω zz p = ω zz m . Closer to the transition, intensity between ω = 0 and the spin-wave peak grows; therefore the characteristic frequency ω xx m becomes smaller than the spin-wave frequency ω p . Their difference is thus a measure of the relative weight of non-single-spin-wave contributions to S(q, ω). Figure 6(a) shows the situation at T KT , where S xx exhibits large non-single-spin-wave contributions and ω xx m < ω xx p . Above the transition, ω xx m is no longer linear in momentum for small q, as shown in fig.  6(b), and differs strongly from the zz-component. The latter still has both a spin-wave peak that is linear in momentum, as well as intensity at small ω, so that ω zz m is smaller than ω zz p . The central question of critical dynamics is that of scaling, i.e. whether data from lattices of different size match when properly scaled. As mentioned in section (2.2), we can test scaling and extract the dynamic critical exponent z in two ways, by analyzing the characteristic frequency ω m , or by looking at S(q, ω) itself.
We concentrate on the dynamic critical behavior of our model at T ≤ T KT , the critical region, in which the correlation length in an infinite system is divergent. The relevant length scale on a finite lattice is therefore the lattice size L, and we expect scaling for suitable functions of qL, as described in section (2.2). From the analytical results we expect the dynamic critical exponent to be z = 1 (eq. (19)).
In figure 7, (a),(b) we show ω xx m L z as a function of qL, for z = 1.00 and at temperatures T ≤ T KT . From eq. (13) we expect the data to fall on a single curve if scaling holds. This is indeed exactly the observed behavior at all temperatures T ≤ T KT . The asymptotic behavior for large L is strictly linear, ω m L z ∼ qL; i.e. for z = 1, ω m ∼ q. each finite lattice size the dispersion curve flattens when q becomes large. Therefore as L increases, the data in figure 7 start to move away from the asymptotic behavior at progressively larger values of qL. We have also analyzed the scaling behavior for different values of z, e.g. z = 1.10 (not shown in the figures). In that case the data for different L diverge from each other immediately; they do not fall onto a common line even at the smallest momenta.
The scaling curves for ω xx m at all three temperatures T ≤ T KT are very similar, with variation only in their slope. In contrast, we do not observe similar scaling behavior in ω xx m at T = 0.8 above the transition (not shown in the figures).
Analyzing the out-of-plane characteristic frequency ω zz m , we found that (for q = 0) it has the same scaling behavior as the in-plane component. At T = 0.4 the data for ω zz m and ω xx m are indistinguishable. When intensity below the spin-wave peak grows in S xx at larger T , the scaling curve for ω xx m has a smaller slope than ω zz m . Interestingly, at T = 0.8, above the transition, not only are there spin-wave peaks present in S zz , but ω zz m also shows the same scaling behavior as below the transition, with z = 1.0.

Finite size scaling of S(q, ω)
If dynamic finite size scaling holds, then the scaled neutron scattering function itself should fall onto a single curve for sufficiently large lattices. Corresponding to eq. (12), figure 8 shows S xx (q, ω)/(L z S(q)) versus ω L z .
We see that at all temperatures T < ∼ T KT the data do indeed fall onto a single line within error bars, when we choose z = 1.00. This is not only true for the spin-wave peak itself, but for the whole range of ω L z . Only the data from very small lattices (not shown here) deviate systematically. Even at T = 0.725 the data scale quite well for the values of L for which we have data. The correlation length at T = 0.725 is still very large (see appendix); deviations from scaling could presumably be seen if data for much larger lattices were available.
Note that scaling with ω L z implies that at fixed qL and for large lattices the spinwave peak is very narrow in units of ω. Its width is therefore very sensitive to the time cutoff in the spin dynamics integration, and we had to use the very long time integrations described in section 3 in order to obtain asymptotic results.
The finite size scaling behavior is very sensitive to variations in z. As an example, figure 8(e) shows that at T = 0.4 the data do not scale when choosing z = 1.05, even upon excluding all lattice sizes L < 96. Using similar plots, we obtain z = 0.99 (2) at T = 0.4 , z = 1.00 (2) at T = 0.6 , z = 1.00 (6) at T = 0.7 , z = 1.02 (3) at T = 0.725 .
(The relatively large error for T = 0.7 is a consequence of the limited amount of data available at this temperature.) It is remarkable that the dynamic critical exponent is the same accross this range of temperatures, whereas the static exponent η varies strongly, from η = 0.082(2) at T = 0.4 to η = 0.247(6) at T = 0.7 (see appendix).
The zz-component of S(q, ω) is extremely narrow at T = 0.4 and T = 0.6, and cannot show scaling given our maximum integration times. At T = 0.7, the spin-wave peak in S zz has become wider, and we do observe scaling, as shown in figure 8(f).
Above the phase transition, the relevant length scale is the correlation length ξ, not the lattice size L. We saw earlier ( fig. 2(d),(e)) that finite size effects are already small for our lattice sizes. Yet, surprisingly, there is scaling-like behavior for small momenta in S xx even above the transition, as show in fig. 9(a). Note that, at constant n q , the horizontal scale ωL is proportional to ω/q. The data do not scale when different n q are compared. At large momenta, a spin-wave peak is visible (see also fig. 3(f)). For the out-of-plane component S zz ( fig. 9(b)), the spin-wave heights do not obey the scaling equation, eq. (12).

Tests of Nelson-Fisher scaling form; large ω-behavior
Nelson and Fisher [8] predicted the scaling form eq. (16) for S xx (q, ω) : This provides an explicit opportunity to compare data at different temperatures and at different values of n q . Equation (16) implies the finite size scaling equation (12) used in the previous section, with z = 1: It also implies for large enough lattice sizes L, for which eq. (16) is valid. (Note that when n q = qL 2π is constant, the arguments ωL in eq. (29) and ω cq in eq. (30) are equivalent.) In fig. 10(a) we use eq. (30) to compare the in-plane scattering function S xx at different temperatures, for constant n q = 2. The same data appear unscaled in fig. 1(a), and with dynamic finite size scaling in fig. 8. Obviously, eq. (30) is not satisfied: the scattering function at different temperatures within the KT-phase does not scale to the same shape Ψ( ω cq ). This is also the case at other values of n q . Different values of n q are compared in fig. 10(b), at T = 0.7. Again, the data do not scale. Moreover, the dependence of the spin-wave peak on n q is not monotonous: for increasing n q the peak height first grows, is approximately constant for n q = 3...8, and then shrinks. As n q becomes large, there is growing intensity below the spin-wave peak (see also fig. 5(d)). Note that equations (16) and (30) are at odds with the fact that the additional peaks we observed in S xx (section 4.2) have positions which do depend on n q . The data in fig. 10  For large frequencies, S(q, ω) appears to be independent of lattice size. We show S xx on a log-log scale in fig. 11, with data scaled similarly as in fig. 10. Nelson and Fisher [8] predict that S xx is also independent of momentum and follows a power law ω −ρ , with ρ = 3 − η (eq. (18)). We see somewhat different behavior, which is determined mainly by n q . For n q = 2, fig. 11(a), the data can be fitted with ρ = 3.0(1) at all temperatures. Note the sizeable structure at low T . For larger n q , fig. 11(b), ρ increases. There is also noticeable curvature in ω, with larger ρ at higher ω; ρ also increases slightly with temperature. The exponents in fig. 11(b) range from 3.7(1) at n q = 4 to 5.4(2) at n q = 24. The out-of-plane correlations S zz (not shown) also decay with a power law with momentum-dependent exponents.
Equation (16) also implies In fig. 12 we use eq. (31) to compare data for different momenta q, at the KT phase transition temperature, using η = 0.25 and constant n q . Here the data do scale. This scaling is also implied by dynamic finite size scaling, fig. 8(c), together with a functional dependence S xx (q) ∼ q η−2 . Note that for constant L (instead of constant n q ), we obtain a non-scaling picture indistinguishable from fig. 10(b).

Lineshapes
In figure 13 we compare our results with theoretical predictions for the shape of S xx (q, ω).
We show data at T = 0.7, for L = 128 and q = π/32 , normalized according to eq. (30), and we compare with predictions using η = 1/4, similarly normalized. peak, we tried to find a functional form which does fit these data reasonably well and should thus be an approximation to the actual form, so that our data can more easily be compared to the results of future theoretical calculations. We found that a modified form of eq. (30) works well, namely a widening of eq. (30) with a Gaussian resolution function eq. (27), with δω a free parameter. Around the spin-wave peak we obtained fairly good agreement with our (unconvoluted) data, as shown by the thick dashed line in Fig. 13, which uses δω = 0.01. However, different values of δω in the modified function are necessary to describe the data at different n q . The modified function can of course not describe the additional structure in S xx (q, ω), including the central peak.
The large frequency behavior of the data is shown in the inset in figure 13. The prediction by Nelson and Fisher agrees with the data at large frequencies qualitatively. However, as described in section 4.5, a fit to S xx (q, ω) for large ω results in a different power law exponent than predicted. Not surprisingly, at large frequencies the prediction by Villain, intended for the spin-wave peak divergence, does not describe the data correctly.
We conclude that below and above the transition, the actual lineshape is quite different from the predicted forms, with a much wider spin-wave peak, and a lot of additional structure.
Above the transition, the theoretical predictions eq. (22) and eq. (23) do not describe the data well either. For the in-plane component S xx we see two different regimes in ω. At small ω, it is compatible with a Lorentzian-like peak ∼ (ω 2 + a) −b , but with an exponent b that is momentum-dependent (e.g. b ≈ 1.1(1) at q = π/48, b ≈ 0.43(2) at q = π/16). At large ω, S xx decays with a power law ∼ ω −c , with c = 3.2(2) for small momenta (see fig. 11(a)). The out-of-plane component S zz does not show the predicted central peak at all; instead it exhibits a spin-wave peak.

Conclusions
We have performed the first high precision study of the dynamic critical behavior of the XY -model, at five different temperatures below, at, and above T KT , on square lattices of size up to 192×192. We have determined the critical temperature to be T = 0.700(5) J/k B . Starting from about 1000 equilibrium configurations generated by an efficient Monte Carlo procedure at each temperature and lattice size, we have integrated the equations of motion of the spins to very large times, t max = 400J −1 , and measured space-displaced, timedisplaced correlation functions to compute the neutron scattering function S(q, ω).
At temperatures up to T KT , S(q, ω) exhibits very strong and sharp spin-wave peaks in the in-plane-component S xx . As T increases, they widen slightly and move to lower ω, but remain pronounced even just above T KT . For increasing momentum they broaden and rapidly lose intensity. Well above T KT , the spin-wave peak disappears in S xx , as expected, and we observe a large central peak instead.
In addition to the spin-wave peak, the in-plane component S xx exhibits a rich structure of small intensity, which we interpret to come from two-spin-wave processes. Furthermore, S xx shows a clear central peak, even below T KT , which becomes very pronounced towards the critical temperature. Neither this strong central peak nor the additional structure are predicted by existing analytical calculations.
The out-of-plane component S zz is much weaker than S xx , except for large q. It displays a sharp spin-wave peak at all temperatures, even above T KT . The peak widens with increasing temperature, and only at low T is it consistent with a delta-function shape.
We measure the dispersion relation, i.e. the position of the spin-wave peak as a function of momentum, to be linear to high accuracy. Its slope, the spin-wave velocity, decreases with increasing temperature approximately linearly, as expected from approximate analytical results at small T.
Examining dynamic finite size scaling, we show that both the characteristic frequency ω m and the neutron scattering function S(q, ω) itself scale very well for all T ≤ T KT , with a dynamic critical exponent of z = 1.00(4) that does not depend on temperature, whereas the static exponent η varies strongly.
The shape of the scaling function is not well described around the spin-wave peak by the available theoretical predictions, nor is the shape of the scattering function above the transition, and the additional structure had not been predicted at all. The data which we have presented here are of sufficiently high quality that meaningful comparison with theory and experiment is possible. We hope that this spin dynamics study will thus serve to stimulate further effort in this area.

Appendix: Static critical behavior
The determination of the transition temperature in XY -like systems has been notoriously difficult [2,3]. The best previous estimate [2] for T KT for the model considered here was T KT = 0.725 ± 0.010, estimated from the onset of vortex-pair unbinding, which is a procedure that is quite difficult to apply with high precision. The results of our high resolution spin dynamics study at T = 0.725 prompted us to perform a new, more accurate determination of T KT , using the powerful hybrid Monte Carlo algorithm described in section 3.
We carried out a set of static Monte Carlo simulations, with lattice sizes L = 64, 128, and 256, and 40000 hybrid Monte Carlo sweeps in each run. We then analyzed the static correlations C(r) = 1 2 S x (0)S x (r) + S y (0)S y (r) in three different ways: (i) using finite size scaling, (ii) with a fit to a power law decay, and (iii) with a fit to the free lattice propagator. The results of all three methods are in excellent agreement.
The finite size scaling ansatz for the correlation function is [28] C(r, L) = r −η f (tL 1/ν , r L ) , wheret is the reduced temperature (t = |1−T /T c |), and ν the correlation length exponent. Since our model is critical throughout the KT-phase, we havet = 0 for all T ≤ T KT . With the correct value of η, the data for different lattice sizes should therefore coalesce on a plot of C(r) L η versus r L . Figure 14 shows such plots at T = 0.7, for η = 0.24, 0.25, and 0.26. Choosing η = 0.25 results in very good scaling over the whole range of distances r. At lower temperatures, scaling of similar quality is achieved with smaller values of η. At a slightly higher temperature of T = 0.71, however, we observed only mediocre scaling, with an effective η of 0.29(1), but with small systematic deviations from scaling already visible. At T = 0.725, the deviations from scaling are stronger still.
Independent estimates of η were obtained from simple power law fits C(r, L) ∼ r −η , for lattice size L = 256, and r max ≤ 20. The results for η depend only very little on the fit-range, and agree within error bars with those from finite size scaling.
A simple power law ansatz ignores the fact that C(r, L) is actually periodic in r with period L. The full functional form can be deduced from the fact that within the KTphase, the model is thought [1,29] to behave like a free field theory, for which the exact finite lattice propagator (correlation function) at an effective "temperature" 1/η is [29] C(r, L) = exp(−η Γ(r, L)), with . (34) We also used this functional form to fit C(r, L). The results for different lattice sizes and different fit ranges (up to r = L/2 and excluding r < 5) agreed with each other and with the results from finite size scaling. The quality of fits was very good for all T ≤ 0.70, whereas for larger temperatures (T = 0.71 and above) it deteriorated strongly, and the results became lattice size dependent. We have found no evidence for logarithmic corrections. It is possible that small corrections are present and introduce a subtle bias. In this case our results for η would need modification. Assuming that logarithmic corrections are indeed negligible, we obtain as our combined results from all three methods η = 0.082 (2) at T = 0.4 η = 0.153 (5) at T = 0.6 η = 0.247 (6) at T = 0.7 .
We note that the expected linear dependence of η on T below T KT does not seem to be satisfied at these temperatures in our model. Both from assuming η = 1/4 at T KT [1], and from the different qualitative behavior of C(r, L) for T ≥ 0.710, we conclude that the Kosterlitz Thouless transition temperature is T KT = 0.700(5)J/k B .
This estimate is slightly below the value from [2] which we had used at the start of the spin dynamics study, and clearly below the estimate of ref. [30]. Slightly above the KT-transition, at T = 0.725, the dynamic behavior of S(q, ω) in our study resembles that at T KT (whereas at T = 0.8 it is very different). This may be explained by looking at the correlation lengths: The correlation function at T = 0.725 can be fitted by an Ornstein-Zernike form Γ(r) ∼ r − 1 2 e −r/ξ with a value of ξ = O(400). Since this correlation length is larger than the lattice sizes we have used in our study, a behavior resembling the KT-phase (where ξ = ∞) is not surprising. A similar (approximate) fit at T = 0.8, on the other hand, gives ξ = O(10), much smaller than our lattice sizes. After our work was completed, we received a paper by Cuccoli et al. [31] with a Monte Carlo study of the statics of our model. Their results are in excellent agreement with ours.