Short-range correlations in modified planar rotator model

We introduce a model inspired from statistical physics that is shown to display flexible short-range spatial correlations which are potentially useful in geostatistical modeling. In particular, we consider a suitably modified planar rotator or XY model, traditionally used for modeling continuous spin systems in magnetism, and we demonstrate that it can capture spatial correlations typically present in geostatistical data. The empirical study of the spin configurations produced by Monte Carlo simulations at various temperatures and stages in the nonequilibrium regime shows that their spatial variability can be modeled by the flexible class of Mat\'{e}rn covariance functions. The correlation range and the smoothness of these functions vary significantly in the parameter space that consists of the temperature and the simulation time. We briefly discuss the potential of the model for efficient and automatic prediction of spatial data with short-range correlations, such as commonly encountered in geophysical and environmental applications.


Introduction
Modern remote sensing techniques allow recording enormous amounts of spatial data, including geographical, natural resources, land use, and environmental remote sensing images, which need efficient (preferably real-time) processing [1,2]. Such processing involves the reconstruction of missing data that often occur due to different reasons, such as equipment malfunctions and gaps in the coverage of the targeted area that appear as a result of restricted satellite paths or bad weather conditions [3,4,5,6,7]. However, such massive data sets can not be efficiently handled by standard geostatistical methods, such as kriging [8]. The main drawbacks of such methods are high computational complexity, difficulties to automatize the algorithms to work without subjective user inputs (selection of variogram and search radius for interpolation) and often the necessity of some data pre-processing (e.g., lognormal transformation) if the data does not comply with the Gaussianity requirement [9]. Thus, there is a need to develop new spatial prediction techniques that overcome these shortcomings [10,11]. Recently, we have proposed efficient spatial classification methods based on models inspired from statistical physics, in particular discrete spin Ising, Potts, and clock models, employing a heuristic "energy matching" principle [12,13].
In the present study we extend the idea of using spin models for spatial prediction purposes and introduce a new method that is based on a continuous spin planar rotator model. Even though the above mentioned nonparametric discrete-spin-based classification models have been shown to be rather competitive, their performance was based on the assumption of the existence of relevant correlations at some unknown parameters and the prediction results were discrete values even for continuous data. As we empirically demonstrate, the modified planar rotator model inherently displays flexible short-range spatial correlations that vary significantly over the model's parameter space and could be used for spatial interpolation of continuous geostatistical data.

Model and simulations 2.1. Standard planar rotator model
The Hamiltonian of the standard two-dimensional planar rotator model with nearest-neighbor interactions on a square lattice is defined as where s i = (cos φ i , sin φ i ) is a continuous spin on i-th lattice site, represented by a twodimensional unit vector, φ i ∈ [0, 2π] is an angle associated with the spin s i , J is an exchange interaction parameter and i, j denotes the sum over nearest neighbors. The Mermin-Wagner theorem [14] prevents any long-range ordering at finite temperatures in such a model, which has, however, been intensively studied in connection with quasi long-range ordering (the so called Kosterlitz-Thouless phase) that appears at low temperatures [15,16]. The phase transition is produced by the unbinding of vortex-antivortex pairs at the Kosterlitz-Thouless critical temperature T KT , below which all spins are almost aligned even though true long-range order is destroyed by spin fluctuations. The quasi long-range ordering for T < T KT is characterized by the power-law decaying correlation function where h is the spin separation distance, L is the linear lattice size, and η(T ) = T /2πJ is the temperature-dependent exponent.

Modified planar rotator model
Even though the algebraic correlations with the tunable exponent η(T ) could be relevant for modeling of some spatial processes, there are some deficiencies that prevent the straightforward use of the standard planar rotator model for such purposes. First of all, one needs to define a proper mapping between the spin values and the geostatistical values. Geostatistical data are often positively correlated reflecting a degree of spatial continuity, which means that neighbors with similar values are more likely (i.e., have lower energy) than those with different values. Apparently, this condition is not fulfilled in the standard planar rotator model, described by the Hamiltonian (1); the latter, for example, assigns equal energies to a pair of neighbors whose angles differ by θ as to a pair whose angles differ by 2π + θ. This degeneracy could be fixed by defining the energy functional so that it monotonically increases with the turn angle between the neighboring spins on the entire interval [0, 2π]. The latter can be achieved by modifying the Hamiltonian (1) to take the following form of the modified planar rotator (MPR) model where 0 < q ≤ 1/2 is a modification factor. In the following we will use a fixed value of q = 1/2. The difference between the nature of spatial correlations in the standard and modified models is apparent from snapshots of spin configurations (turn angles) in the low-temperature region, shown in Fig. 1. While the presence of the vortex-antivortex pairs shows up in Fig. 1(a) in the form of sharp boundaries between the domains of similarly oriented spins, the variation of the spin values in the MPR model, shown in Fig. 1(b) is rather smooth, as it is typical for geostatistical data.

Monte Carlo simulations
We use Monte Carlo (MC) simulations with the Metropolis update rule and vectorized checkerboard algorithm. The results are presented for a square lattice with the size L × L.
We chose L = 128 as a compromise value between the computational resources needed for MC simulations and calculations of the empirical variograms (see below) on one side and the effort to secure ergodic conditions on the other side. In typical simulations of magnetic systems one would like to suppress boundary effects by employing periodic boundary conditions (see the snapshots in Fig. 1). However, these are not appropriate for simulation of geostatistical data, since normally there is no reason to assume the data on the opposite boundaries are correlated. Free boundary conditions are not appropriate either since the spatial process is not necessarily confined within the domain boundaries. Instead, we consider it reasonable to apply boundary conditions that assume the smooth continuation of the spatial process beyond the borders. This assumption is implemented by inserting additional nodes that decorate the lattice and requiring that these additional points have the same values as their nearest neighbors inside the lattice. Therefore, if s i,j is a spin in the i-th row and j-th column of the lattice, i, j = 1, . . . L, then 1 and s 0,j = s 1,j , where the rows and columns with indices 0 and L + 1 respectively are formed by auxiliary nodes around the lattice added to deal with the boundary effects.

Modeling spatial variability
The MC simulations produce spatially correlated spin realizations s i , which can be represented by their turn angles φ i . It turns out that there is great flexibility in the spatial correlations as we move in the temperature and simulation time parameter space. Therefore, in order to model the spatial variability (the local variogram) of the simulated spin values we use a very flexible Matérn covariance model. In addition to the trivial parameters of the standard models (Gaussian, exponential, spherical, etc.), controlling the variance, σ 2 , and the characteristic covariance length, ξ, the Matérn model involves one more parameter, ν, controlling the smoothness of the spatial process. Therefore, it can be considered appropriate for those geostatistical applications in which controlling the smoothness is important. Furthermore, the Matérn model includes the Gaussian and the exponential models as special cases for ν → ∞ and ν = 0.5, respectively, as well as several other bounded models [17]. Due to its flexibility, it has been used in various areas including pedology [17,18,19], hydrology [20,21], topography [22], health modeling [23], meteorology [24] and environmental modeling [25]. The Matérn correlation function has the general form where θ ′ = (ξ, ν) and K ν is the modified Bessel function of order ν. Then the corresponding variogram function γ(h; θ), which is typically used in geostatistics, is related to the correlation function as where σ 2 is the random field variance, σ 2 n is the nugget variance that corresponds to uncorrelated fluctuations and θ = (σ n , σ, ξ, ν) is a vector of the complete model parameters.
Having generated realizations of spin values from MC simulations, we can assess their spatial variability by calculating the experimental (i.e., sample-based) variogram aŝ where N (h) is the number of pairs of spins separated by the distance h.
There are several methods of fitting the model to the experimental variogram. In the present study we used the weighted least squares estimator (hereafter WLS) that was reported to give the best overall results [27]. The WLS estimator is based on minimizing the weighted sum of squared residuals (objective function), which is defined as follows: whereγ(h k ) is the experimental and γ(h k ; θ) the model variogram. The summation runs over the k = 1, . . . , k max lags containing N (h k ) pairs of points, where h(k max ) is less than one half of the maximum lag in the sample.

Results and discussion
In order to empirically study spatial correlations in the realizations of spin values generated by the modified planar rotator model we run MC simulations at various temperatures T (in units J/k B , where k B is Boltzmann constant). Each simulation starts from an initial configuration of spatially uncorrelated random spin angles uniformly distributed in the interval [0, 2π] and the snapshots are collected after τ MC sweeps. Then, using Eq. (6) the experimental variogram is calculated and subsequently fitted to the appropriate model by the WLS method defined by Eq. (7). The results obtained at relatively higher temperatures are demonstrated in Fig. 2, for T = 10 −2 and τ = 10, 100, 1000 and 10000. Soon after the start of the simulation one can observe the development of spatial correlations both in the snapshots in the form of nucleation of small spin domains with similar values as well as in the small-scale behavior of the experimental variograms. The fitted parametersξ = 1.62 andν = 0.62 for τ = 10 in Fig. 2(a) indicate that the realization is quite rough with a rather small characteristic length. As the relaxation proceeds the variance decreases and both the characteristic length and the smoothness parameter increase, as shown in Figs. 2(b) and 2(c), for τ = 100 and 1000, respectively. Nevertheless, as equilibrium is approached the realizations become rougher again, whereas the characteristic length further increases and saturates only in equilibrium to a temperature-dependent value. The estimated characteristic lengthξ = 56.51 for the equilibrium configuration at τ = 10000 in Fig. 2(d) implies that the samples on the 128 × 128 lattice suffer from lack of ergodicity, which is also reflected in the trend displayed by the corresponding experimental variogram.
As the temperature is decreased the trend of building up correlations in the equilibration process is further enhanced. Since thermal fluctuations are gradually suppressed the realizations become smoother. This tendency is apparent by visual comparison of the snapshots obtained at  Figure 3. The same as in Fig. 2 for T = 10 −4 . Expo(θ) and Gaus(θ) denote respectively the exponential and Gaussian models for the inferred parameter valuesθ = (σ n ,σ,ξ). Note that in Fig. 3(a) the fitted Matérn and exponential models coincide as do in Fig. 3(d) the Matérn and the Gaussian models. T = 10 −4 , shown in Fig. 3 with those generated in the same stages of the relaxation at T = 10 −2 , shown in Fig. 2. The difference between the inferred values of the smoothness parameter for the respective cases becomes evident particularly at later stages (τ = 10 3 and 10 4 ). As mentioned above, the exponential covariance model with quite rough (non-differentiable) realizations and the Gaussian model with very smooth (infinitely differentiable) realizations can be obtained from the Matérn model as special cases for ν = 0.5 and ν → ∞, respectively. To emphasize the flexibility of the correlations developed in the relaxation process we also fit the relatively rough data for τ = 10 ( Fig. 3(a)) to the exponential and the very smooth data for τ = 10 4 (Fig. 3(d) (see the explanation below), the characteristic length is expected to increase with decreasing temperature and increasing simulation time. This expectation can be justified by the fact that the ground state of the MPR model (the equilibrium state at T = 0 corresponding to the lowest energy), is the ferromagnetic state with all the spins aligned in the same direction. Thus, for T → 0 one can expect that the correlation (characteristic) length will span the entire lattice and eventually go to infinity for L → ∞. Next, we comment on the big discrepancy between the characteristic length parametersξ, estimated for the Matérn and Gaussian models, for T = 10 −4 and τ = 10 4 . In fact, in the Matérn model the parameter ξ has been found to be highly negatively correlated with the smoothness parameter ν [18,26,28]. This correlation is also apparent by looking at the objective function of the WLS fit, defined by Eq. (7), which is shown in Fig. 4 in the ξ−ν parameter space, for the data presented in Fig. 3(d). It demonstrates that the samples with smallξ and largeν values are quite "close" to the samples with largeξ and smallν values. Put differently, a small difference between samples from the same population (particularly in non-ergodic samples) can lead to completely different parameter estimates. Thus, one should be careful when performing parameter inference from the empirical variogram. Recently, we have presented a model-independent definition of the correlation length borrowed from statistical field theory and proposed to use a so-called ergodicity index to compare coarse-grained measures corresponding to both trivial (standard) and non-trivial, e.g. Matérn, covariance models with different parameters [28].

Summary and outlook
We have modified the standard planar rotator spin model from statistical physics to display a flexible type of short-range spatial correlations relevant in geostatistical modeling. In particular, the empirical study of spin configurations produced by Monte Carlo simulations in the nonequilibrium regime at various temperatures and stages shows that the smoothness and correlation range vary greatly versus the temperature and the simulation time. This behavior implies that the model has good potential for the simulation and prediction of spatial processes in geophysical and environmental applications. In particular, one can use the model to perform conditional simulations of gridded data, such as remote sensing images. Conditional simulations honor the sample values and reconstruct the variability of missing data, based on simulation parameters inferred from the available samples. Owing to the fact that the model does not show undesirable critical slowing down, the relaxation process is rather fast; furthermore, the short-range nature of the interactions between spin variables allows vectorization of the algorithm. Consequently, the proposed method is significantly more efficient than the conventional geostatistical approaches, and thus applicable to huge datasets, such as satellite and radar images. The implementation details are now under investigation.