This is used to compute the dynamical structure factor, also called scattering function:
\begin{equation*} S(\vec Q, t) = \langle s-\vec Q(0) ⋅ s\vec Q(t) \rangle, \end{equation*} with \begin{equation*} s\vec Q(t) = ∑i, α si, α(t) e^{-i (\vec R_i + \vec r_α) ⋅ \vec Q} \end{equation*}
The problem is twofold. First, it is a dynamical problem, that is it requires to solve the equations of motion of the system to see how each spin evolve with time. Second, it is a statistical problem, in the sense that the scattering function is defined as an average over the canonical ensemble (fixed temperature \(T\)). This is why it will be solve in a two steps process:
- Sample from the canonical ensemble using MC methods
- Evolve each sample in time using semi-classical methods
The single spin-flip Metropolis Hastings algorithm is used. However, it can be quite inefficient at low temperature due to the number of rejected samples. To reduce this effect, the solid angle for each spin flip is reduced to keep the acceptance rate above 0.4. A pseudocode implementation is given below:
Let S[i] the initial system state. Repeat N_steps times: While not accepted: Choose a random site i. Choose a random 3D unit vector S'. Let E the energy of the state S[i]. Let E' the energy of the state S[i], with S[i] replaced with S'. If E' < E: Accept. Else: Accept with probability exp((E - E') / kT) End If End While End Repeat
How is actually implemented the solid angle constrained practically ? Having a running estimator of the acceptance rate, and updating the solid angle accordingly ?
- Number of samples: 1000
- Stride between two samples: enough such that the stochastic correlation is less than 0.1.
how is actually defined the stochastic correlation ? \(Corr(\{\vec S_i\}, \{\vec S’_i\})\) ?
At low temperature, the system has trouble to get out of the local minimum of energy to explore other regions of the phase space. To mitigate this issue, an overrelaxation scheme is used. Recall that the hamiltonian can be written in terms of the local field \(h_α(\vec R)\) at each site \((\vec R, α)\): \begin{equation*} H = ∑\vec R \vec R’, α, β Jα, β(\vec R’ - \vec R) \vec S_α(\vec R) ⋅ \vec S_β(\vec R’) = ∑\vec R, α \vec h_α(\vec R) ⋅ \vec S_α(\vec R) \end{equation*} One can see that the energy does not change if a spin is rotated, while keeping the angle with its local field fixed. This is exploited by the overrelaxation method. Each time a spin is selected, it is first rotated around its local field, then assigned a random value. If it is rejected, the rotation is still kept.
Is it correct ? Especially, is it applied to each selected spin, or only to failing ones ?
Using the Heisenberg picture equations of motion \(\pdiff{A}{t} = i/\hbar [H, A]\), and the spin 1/2 commutation relation \([\hv S^i, \hv S^j] = εijk \hv S^k\), one can obtain the (semiclassical) non-linear Bloch equations: \begin{equation*} \diff{s_i(t)}{t} = - s_i(t) × \left( ∑_j Jij s_j(t) \right) \end{equation*} This equation is simply numerically integrated. An 8th order explicit Runge-Kutta scheme is used.
The RK error parameter as well as the RK order have been fixed in order to preserve the Euclidean distance \(d = [∑_i (\vec s_iRK - s_jBS)^2]1/2\), i.e. the distance between time trajectories obtained with the RK method and with the more robust but time consuming Burlisch-Stoer (BS) algorithm.
I could not find how the \(\frac{8(8-1)}2 = 28\) RK8 params are defined, along with the timestep.