MARQOV is a yet to be released software framework, currently developed at the Institute for Theoretical Physics and Astrophysics at the University of Würzburg. The acronym stands for MAssively paRallel Quenched disOrder on Variable geometries. The particular focus on efficient and highly scalable Markov-Chain Monte Carlo simulations of classical equilibrium spin models on frustrated and disordered geometries renders MARQOV unique in the world of open-source research software.

The General Hamiltonian

The general Hamiltonian builds the centerpiece of MARQOV, and allows to cover a large number of physical models and general geometries. It governs the physical interactions in the system and is formally split into three parts

\[ \mathcal{H} = \mathcal{H}_\text{int} + \mathcal{H}_\text{self} + \mathcal{V}_\text{pot}, \]

interaction terms, on-site single particle contributions and contributions due to a potential, which read

\[ \begin{aligned} \mathcal{H}_\text{int} &= \sum\limits_{\alpha=1}^{N_\alpha}J^{(\alpha)} \sum\limits_{\alpha\text{-bonds: }i,j} \phi_i \, D^{(\alpha)}_{ij}(\phi_i,\phi_j)\, \phi_j\\ \mathcal{H}_\text{self} &= \sum\limits_{\beta=1}^{N_\beta}h^{(\beta)} \sum\limits_{\beta\text{-sites: }i} g^{(\beta)} (\phi_i)\\ \mathcal{V}_\text{pot} &= \sum\limits_{\gamma=1}^{N_\gamma}k^{(\gamma)} \sum\limits_{\gamma\text{-sites: }i} v^{(\gamma)} \big(\phi_i,\ldots\big) \end{aligned} \]

The configuration (or state) of the individual spins, which are placed on discrete lattice positions \(i\), is denoted by the symbol \( \phi\)1. This covers, e.g. Ising spins where \(\phi\in\pm 1\) and \(N\)-vector models, where \(\phi\in\mathbb{R}^N\). The set of all individual spin configurations hence constitutes the state space of the physical system.

Note that the first term in the above Hamiltonian does only support interactions which can be cast into the given matrix-vector form. More general two-site interactions \(f(\phi_i,\phi_j)\) are not natively covered. However, we provide convenient interfaces where users can nevertheless implement more general interaction terms, though, at the cost of not being able to access the optimized set of standard Monte Carlo update algorithms that come along with MARQOV.

The on-site term allows for full flexibility in implementing self-interactions of any functional form, with the only requirement that the dimension of the target set of global pre-factor \(g^{(\beta)}\) needs to match the dimension of the global constant \(h^{(\beta)}\) in order to provide a meaningful scalar product. Compare the XXZ Heisenberg model below for an example.

The potential term is used to account for any physical interactions which go beyond two-site and self-interactions. Typical examples are long-range Coulomb interactions in a spin ice system or plaquette interactions.

One particular strength of MARQOV is that it is designed for supporting irregular geometries. This is introduced to the above Hamiltonian by the Greek letters, indicating distinct families of bonds \(\alpha\) and spins \(\beta\), \(\gamma\), respectively. Let us clarify this concept using the following figure. In panel (a) there is only one type (family) of sites (open circles), but two families of bonds, in this case referring to nearest (black lines) and next-nearest neighbor connections (blue lines). In the second panel (b) we see three different species of sites (red, blue, green), translating to three \(\beta\)-families in our framework. Finally, panel (c) shows a typical example of a disordered lattice, with only one family of either bonds and sites.



We present a selection of Hamiltonians which fit into this framework with the respective substitutions given:

  • Example 1: Ferromagnetic Ising model in homogeneous external field

    In the celebrated Ising model, individual spins are binary variables. It represents the most basic toy model for ferromagnetic behavior if the interaction constant \(J\) is set positive.

    \[ \mathcal{H} = -J\sum\limits_{\langle i,j\rangle} \phi_i\phi_j + h\sum\limits_{i} \phi_i \]

    where \(\phi_i=\pm 1\). In our framework this translates to one interaction term and one on-site term, which exhibit rather simple structure:

    \[ N_\alpha=N_\beta = 1, \quad J^{(1)} = -J, \quad h^{(1)} = h \]

    \[ D_{ij} = \mathrm{id}, \quad g^{(1)} = \phi_i \]

  • Example 2: Diluted Heisenberg model

    In this model, spins are three-dimensional unit vectors, placed on a regular lattice with quenched random vacancies or impurities. The system obeys a global \(O(3)\) symmetry, which means that the Hamiltonian is invariant under three-dimensional rotations in spin space. It is formally given by:

    \[ \mathcal{H} = J\sum\limits_{\langle i,j\rangle} \epsilon_i \epsilon_j \vec{\phi_i}\vec{\phi_j} \]

    \[ \text{where}\quad \epsilon_i,\epsilon_j\in\{0,1\} \]

    \[ \text{and}\quad \vec{\phi}_i\in\mathbb{R}^3, \quad |\vec{\phi}_i|^2=1. \]

    Here, binary variables \(\epsilon\) are used to encode whether the corresponding site is magnetic (\(\epsilon=1)\) or a non-interacting impurity (\(\epsilon=0)\). Note that there is some redundancy in this notation. If the lattice, which is the set of sites \(i\) and the set of bonds \((i,j)\), does not include impurity sites in the first place, the \(\epsilon\)'s can be discarded.

    \[ N_\alpha=1, \quad J^{(1)} = J, \quad D_{ij} = \mathrm{id} \]

  • Example 3: Potts model

    The \(q\)-state Potts model represents a generalization of the Ising model, where spins can take on (q values, rather than only two. The Hamiltonian is given by

    \[ \mathcal{H} = J\sum\limits_{\langle i,j\rangle}\delta_{\phi_i,\phi_j} \]

    where \( \phi_i\in{0,1,\ldots,q-1} \) and \( \delta \) represents the Kronecker delta. More details can be found in our Gallery.

  • Example 4: Clock model

    The \(q\)-state Clock model reads

    \[ \mathcal{H} = J\sum\limits_{\langle i,j\rangle}\cos({\phi_i-\phi_j}), \]

    where \(\theta_i = 2\pi k/q\) and \(k \in{0,1,\ldots,q-1}\). In 2D it experiences a BKT transition for \(q\geq 5\), whereas the clock model for \(q=4\) comprises two sets of the Ising model and the 3-state clock model is equivalent to the 3-state Potts model. The clock model for \(q=2\) is simply the Ising model. Furthermore, sending \(q\to\infty\) we arrive at the XY or O(2) model.

    Note that the model matches the form of the interaction in our general Hamiltonian, since by means of trigonometric identities, it can be rewritten as

    \[ \mathcal{H} = J\sum\limits_{\langle i,j\rangle} \vec{\phi_i}\vec{\phi_j} \quad\text{where}\quad \phi_i=(\cos\theta_i, \sin\theta_i) \]

  • Example 5: Spin glasses

    If we generalize an Ising or \(\mathcal{O}(N)\) symmetric model to fluctuating local couplings that can take on random positive and negative values, the system behaves neither ferro- nor antiferromagnetic but rather becomes a prototypical spin glass. The Hamiltonian formally reads

    \[ \mathcal{H} = J\sum\limits_{\langle i,j\rangle} \epsilon_{ij} \phi_i \phi_j \]

    \[ \text{where}\quad \epsilon_i,\epsilon_j\in\left[-1,1\right] \text{(random)} \]

    Note that if the disorder is static (quenched), the interaction strengths \(D_{ij}=\epsilon_{ij}\) are typically read from a table.

  • Example 6: Spin ice with dipolar interactions

    The term spin ice refers to substances, in which the disorder of the magnetic moments at low temperatures behaves qualitatively similar to the proton disorder in water ice. The minimal model is the following:

    \[ \mathcal{H} = - J \sum\limits_{\langle i,j\rangle}\phi_i \phi_j + Da^3 \sum\limits_{i>j}\Bigg[\frac{\phi_i\phi_j}{|\vec{r}_{ij}|^3}- 3 \frac{(\vec{\phi}_i\cdot \vec{r}_{ij})(\vec{\phi}_j\cdot \vec{r}_{ij})}{|\vec{r}_{ij}|^5}\Bigg] \]

    This Hamiltonian is typically defined on a pyrochlore lattice. The first term represents a regular nearest neighbor Ising-type interaction. The second term, however, introduces long-range dipolar interactions and enters via a suitable potential energy term in our framework. Note that for this kind of interactions the sum can be efficiently computed using the Ewald summation technique. Formally, we have

    \[ N_\gamma=1, \quad k^{(1)} = Da^3 \]

    \[ v^{(1)} = v^{(1)} \Big( \phi_i,\{\phi_{j\neq i}\},\{\vec{r}_{ij}\} \Big) = \frac{1}{2}\sum\limits_{j\neq i}\Bigg[\frac{\phi_i\phi_j}{|\vec{r}_{ij}|^3} - 3 \frac{(\vec{\phi}_i\cdot \vec{r}_{ij})(\vec{\phi}_j\cdot \vec{r}_{ij})}{|\vec{r}_{ij}|^5}\Bigg]. \]

  • Example 7: XXZ Heisenberg antiferromagnet in Y external field

    We start with a regular \(\mathcal{O}(N)\) symmetric model in an external field, i.e.

    \[ N_\alpha=N_\beta = 1, \quad \phi_i\in\mathbb{R}^3, \quad |\phi_i|^2=1 \]

    This example demonstrates how direction-dependent couplings can be introduced:

    \[ J^{(1)} = -1, \qquad h^{(1)} = (0, h_y, 0)^T, \qquad D_{ij}=\begin{pmatrix} J_x & 0 & 0 \\ 0 & J_x & 0 \\ 0 & 0 & J_z \end{pmatrix} \]

    \[ g^{(1)}:\phi_i\mapsto\mathbb{R}^3, \qquad g^{(1)}(\phi_i) = \phi_i \]

    A related model, the XXZ Heisenberg antiferromagnet with additional single-ion anisotropy can be found in our Gallery.

  • Example 8: Blume-Capel model

    In the Blume-Capel model, the spins take on integer values \(\phi_i\in{-1,0,+1}\) and couple via an Ising-like interaction. Additionally, we have a single-ion anisotropy term (also called reduced local crystal field).

    \[ \mathcal{H} = -J\sum\limits_{\langle i,j\rangle} \phi_i\phi_j + D\sum\limits_{i} \phi_i^2 \]

    In our framework this translates to

    \[ N_\alpha=N_\beta = 1, \quad h^{(1)} = D, \quad D_{ij}=\mathrm{id}, \quad g^{(1)} = \phi_i^2 \]

  • Example 9: Improved O(N) lattice field theory

    An important model featuring a standard interaction and two on-site terms is the scalar lattice field theory Hamiltonian with additional fourth-order self-interaction term. It is given by

    \[ \mathcal{H} = -\beta\sum\limits_{\langle i,j\rangle} \phi_i \phi_j + m\sum\limits_i \phi_i^2 + \lambda\sum\limits_i (\phi_i^2-1)^2 \]

    In our framework this straightforwardly translates to

    \[ N_\alpha=1, \quad N_\beta=2, \quad \phi_i\in\mathbb{R}^N \]

    \[ \quad J^{(1)} = -\beta, \quad h^{(1)}=m, \quad h^{(2)}=\lambda \]

    \[ g^{(1)} = \phi_i^2, \quad g^{(2)} = (\phi_i^2-1)^2 \]

    Note that the families \(\beta=1\) und \(\beta=2\) may of course contain the same sites.

  • Example 10: Ashkin-Teller model

    The \(N\)-color Ashkin-Teller model is given by

    \[ \mathcal{H} = - J\sum\limits_{\alpha=1}^N\sum\limits_{\langle i,j\rangle} \phi^\alpha_i\phi^\alpha_j - K \sum\limits_{\alpha<\beta}\sum\limits_{\langle i,j\rangle} \phi^\alpha_i\phi^\alpha_j \phi^\beta_i\phi^\beta_j \]

    and features a four-particle interaction term.

    At a first glance, an implementation of this Hamiltonian is not straightforwardly possible based on the above general form of the MARQOV Hamiltonian, as it requires an interaction term that couples spins of different families. However, one can use standard embedding techniques to numerical handle the Ashkin-Teller as \(N\) Ising models, coupled via their interaction strengths.2

  • Example 11: Real projective sigma model

    \[ \mathcal{H} = -\frac{J}{2}\sum\limits_{\langle i,j\rangle} \left(\vec{\phi}_i\vec{\phi}_j\right)^2 \]

    Due to the interaction term contributing quadratically, this model also does not fit the general Hamiltonian in a straightfoward way. However, there is again embedding techniques available, as described in detail in Reference 3.

  1. Other symbols used in the literature to describe configurations of classical spins are \(\sigma\) and \(S\).
  2. Zhu, Q., Wan, X., Narayanan, R., Hoyos, J. A., & Vojta, T. (2015). Emerging criticality in the disordered three-color Ashkin-Teller model. Physical Review B, 91(22), 224201.
  3. Caracciolo, S., Edwards, R. G., Pelissetto, A., & Sokal, A. D. (1993). Wolff-type embedding algorithms for general nonlinear σ-models. Nuclear Physics B, 403(1-2), 475-541.