Restricted Boltzmann Machines for Quantum States with Nonabelian or Anyonic Symmetries

TL;DR

This paper introduces a symmetry-respecting neural network approach using basis sampling in irreducible representations to accurately model quantum states with nonabelian and anyonic symmetries, improving over previous basis-dependent methods.

Contribution

It proposes a basis sampling method in irreducible representations for neural network quantum states, overcoming basis dependence issues in nonabelian symmetry modeling.

Findings

Achieved higher accuracy in ground state calculations of the AFH model.

Successfully targeted excited states and calculated energy gaps.

Generalized the framework to anyonic spin chains.

Abstract

Although artificial neural networks have recently been proven to provide a promising new framework for constructing quantum many-body wave functions, the parameterization of a quantum wavefunction with nonabelian symmetries in terms of a Boltzmann machine inherently leads to biased results due to the basis dependence. We demonstrate that this problem can be overcome by sampling in the basis of irreducible representations instead of spins, for which the corresponding ansatz respects the nonabelian symmetries of the system. We apply our methodology to find the ground states of the one-dimensional antiferromagnetic Heisenberg (AFH) model with spin-half and spin-1 degrees of freedom, and obtain a substantially higher accuracy than when using the $s_z$-basis as input to the neural network. The proposed ansatz can target excited states, which is illustrated by calculating the energy gap of the AFH model. We also generalize the framework to the case of anyonic spin chains.

Equations

$$E_{RBM}(\mathbf{x},\mathbf{h};\mathcal{W})=\sum _i^N\sum _j^M{w}_{ij}{x}_i{h}_j+\sum _i^N{a}_i{x}_i+\sum _j^M{b}_j{h}_j.$$

$$\Psi (\mathbf{x};\mathcal{W})\equiv \sum _{\mathbf{h}}{e}^{-E_{RBM}(\mathbf{x},\mathbf{h};\mathcal{W})},\mathcal{W}\subset \mathbb{C}.$$

$$E(\Psi ;\mathcal{W})=\frac{⟨\Psi ⟩\hat{H}\Psi }{⟨\Psi ⟩\Psi },$$

$${\hat{H}}_{AFH}=\sum _{i=1}^{N-1}{\hat{\mathbf{s}}}_i\cdot {\hat{\mathbf{s}}}_{i+1},$$

$$\begin{array}{rl}|s_1{s}_2;j_2{m}_{j_2}⟩=\sum _{m_{s_1},m_{s_2}}& ⟨s_1{m}_{s_1},s_2{m}_{s_2}⟩j_2{m}_{j_2}\\ & |s_1{m}_{s_1},s_2{m}_{s_2}⟩.\end{array}$$

$$|\Psi ⟩=\sum _{j_k}\Psi (j_1,...,j_{N-1})|j_1,...,j_{N-1};J{M}_J⟩,$$

$${\hat{H}}_{GC}=-\sum _{i=1}^{N-1}{\hat{\mathcal{P}}}_{i,i+1}^1,$$

$$\begin{array}{rl}|s_1{s}_2;j_2{m}_{j_2}⟩& =\sum _{m_{s_1}=-s_1}^{s_1}\sum _{m_{s_2}=-s_2}^{s_2}(-1{)}^{-s_1+s_2-m_{j_2}}\sqrt{2j_2+1}\left(\begin{array}{ccc}{s}_1& s_2& j_2\\ m_{s_1}& m_{s_2}& -m_{j_2}\end{array}\right)|s_1{m}_{s_1},s_2{m}_{s_2}⟩\\ & =\sum _{m_{s_1}=-s_1}^{s_1}\sum _{m_{s_2}=-s_2}^{s_2}⟨s_1{m}_{s_1},s_2{m}_{s_2}⟩s_1{s}_2;j_2{m}_{j_2}|s_1{m}_{s_1},s_2{m}_{s_2}⟩.\end{array}$$

$$\hat{1}=\sum _{m_{s_1}=-s_1}^{s_1}\sum _{m_{s_2}=-s_2}^{s_2}|s_1{m}_{s_1},s_2{m}_{s_2}⟩⟨s_1{m}_{s_1},s_2{m}_{s_2}|.$$

$$\begin{array}{rl}& |s_1...s_N{j}_1...j_{N-1};J{M}_J⟩=\\ & \sum _{\{m_{s_i}\}}\sum _{\{m_{j_i}\}}\left(\prod _{i=1}^N(-1{)}^{-j_{i-1}+s_i-m_{j_i}}\sqrt{2j_i+1}\left(\begin{array}{ccc}{j}_{i-1}& s_i& j_i\\ m_{j_{i-1}}& m_{s_i}& -m_{j_i}\end{array}\right)\right){\delta }_{j_{N}J}{\delta }_{m_{j_N}{M}_J}|s_1{m}_{s_1},s_2{m}_{s_2},...,s_N{m}_{s_N}⟩.\end{array}$$

$$\begin{array}{rl}& ⟨j_0^{\prime }{j}_1^{\prime }...j_{N-1}^{\prime };J{M}_J⟩j_0{j}_1...j_{N-1};J{M}_J\\ & \begin{array}{rl}=\sum _{\{m_{s_i}\}}\sum _{\{m_{j_i}\}}\sum _{\{m_{j_i^{\prime }}\}}& \left(\prod _{i=1}^N(-1{)}^{-j_{i-1}+s_i-m_{j_i}}\sqrt{2j_i+1}\left(\begin{array}{ccc}{j}_{i-1}& s_i& j_i\\ m_{j_{i-1}}& m_{s_i}& -m_{j_i}\end{array}\right)\right)\\ & \left(\prod _{i=1}^N(-1{)}^{-j_{i-1}^{\prime }+s_i-m_{j_i^{\prime }}}\sqrt{2j_i^{\prime }+1}\left(\begin{array}{ccc}{j}_{i-1}^{\prime }& s_i& j_i^{\prime }\\ m_{j_{i-1}^{\prime }}& m_{s_i}& -m_{j_i^{\prime }}\end{array}\right)\right)\end{array}\\ & ={\delta }_{j_1,j_1^{\prime }}{\delta }_{j_2,j_2^{\prime }}...{\delta }_{j_{N-1},j_{N-1}^{\prime }},\end{array}$$

$$\sum _{m_{j_{i-1}}}\sum _{m_{s_i}}(2j_i+1)\left(\begin{array}{ccc}{j}_{i-1}& s_i& j_i\\ m_{j_{i-1}}& m_{s_i}& m_{j_i}\end{array}\right)\left(\begin{array}{ccc}{j}_{i-1}& s_i& j_i^{\prime }\\ m_{j_{i-1}}& m_{s_i}& m_{j_i^{\prime }}\end{array}\right)={\delta }_{j_i^{\prime },j_i}{\delta }_{m_{j_i^{\prime }},m_{j_i}},$$

$$\hat{H}=\sum _{i=1}^{N-1}{\hat{\mathbf{s}}}_i\cdot {\hat{\mathbf{s}}}_{i+1},$$

$$F_{s_{i+1},s_i,j_{i-1},j_{i+1}}^{S,j_i}=(-1{)}^{s_i+s_{i+1}+j_{i-1}+j_{i+1}}\sqrt{(2S+1)(2j_i+1)}\left\{\begin{array}{ccc}{s}_{i+1}& s_i& S\\ j_{i-1}& j_{i+1}& j_i\end{array}\right\}.$$

$$|j_0...j_{i-1}{j}_i{j}_{i+1}...j_{N-1};J{M}_J⟩=\sum _S{F}_{s_{i+1},s_i,j_{i-1},j_{i+1}}^{S,j_i}|j_0...j_{i-1}S{j}_{i+1}...j_{N-1};J{M}_J⟩.$$

$$\begin{array}{rl}& \frac{1}{2}(S^2-s_i^2-s_{i+1}^2)|j_0...j_{i-1}{j}_i{j}_{i+1}...j_{N-1};J{M}_J⟩\\ & =\sum _S\frac{1}{2}(S(S+1)-s_i(s_i+1)-s_{i+1}(s_{i+1}+1))F_{s_{i+1},s_i,j_{i-1},j_{i+1}}^{S,j_i}|j_0...j_{i-1}S{j}_{i+1}...j_{N-1};J{M}_J⟩.\end{array}$$

$$\begin{array}{rl}& ⟨j_0...j_{i-1}{j}_i^{\prime }{j}_{i+1}...j_{N-1};J{M}_J⟩{\hat{\mathbf{s}}}_i\cdot {\hat{\mathbf{s}}}_{i+1}{j}_0...j_{i-1}{j}_i{j}_{i+1}...j_{N-1};J{M}_J\\ & =\sum _S\frac{1}{2}(S(S+1)-s_i(s_i+1)-s_{i+1}(s_{i+1}+1))F_{s_{i+1},s_i,j_{i-1},j_{i+1}}^{S,j_i}{F}_{s_{i+1},s_i,j_{i-1},j_{i+1}}^{S,j_i^{\prime }}.\end{array}$$

$$\hat{U}(\mathbf{\theta })=e^{i({\theta }_x{\hat{s}}_x+{\theta }_y{\hat{s}}_y+{\theta }_z{\hat{s}}_z)}=e^{i\hat{\mathbf{s}}\cdot \mathbf{\theta }},$$

$${\hat{\sigma }}_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\hat{\sigma }}_y=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),{\hat{\sigma }}_z=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).$$

$${\hat{s}}_x=\left(\begin{array}{ccc}0& \frac{1}{\sqrt{2}}& 0\\ \frac{1}{\sqrt{2}}& 0& \frac{1}{\sqrt{2}}\\ 0& \frac{1}{\sqrt{2}}& 0\end{array}\right),{\hat{s}}_y=\left(\begin{array}{ccc}0& \frac{-i}{\sqrt{2}}& 0\\ \frac{i}{\sqrt{2}}& 0& \frac{-i}{\sqrt{2}}\\ 0& \frac{i}{\sqrt{2}}& 0\end{array}\right),{\hat{s}}_z=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right).$$

$$|s_1{s}_2;j{m}_j⟩=\sum _{m_{s_1},m_{s_2}}⟨s_1{m}_{s_1},s_2{m}_{s_2}⟩j{m}_j|s_1{m}_{s_1},s_2{m}_{s_2}⟩,$$

$$e^{i\hat{\mathbf{j}}\cdot \mathbf{\theta }}|s_1{s}_2;j{m}_j⟩=\sum _{m_{s_1},m_{s_2}}⟨s_1{m}_{s_1},s_2{m}_{s_2}⟩j{m}_j{e}^{i({\hat{\mathbf{s}}}_1+{\hat{\mathbf{s}}}_2)\cdot \mathbf{\theta }}|s_1{m}_{s_1},s_2{m}_{s_2}⟩.$$

$$C_{\beta ,\gamma }(i,j)=⟨{\hat{s}}_i^{\beta }{\hat{s}}_j^{\gamma }⟩,$$

$$⟨J^{\prime }{M}_J^{\prime }⟩{\tilde{T}}_{l,m_l}J{M}_J=⟨J^{\prime }{M}_J^{\prime };l{m}_l⟩J{m}_J⟨J^{\prime }⟩{\tilde{T}}_{l}J.$$

$${\tilde{C}}_{0,0}(i,j)=-\frac{1}{\sqrt{3}}\left(C_{x,x}(i,j)+C_{y,y}(i,j)+C_{z,z}(i,j)\right).$$

$$\begin{array}{rl}& C_{x,x}(i,j)=-\frac{1}{\sqrt{3}}\left({\tilde{C}}_{0,0}(i,j)-\frac{1}{\sqrt{6}}{\tilde{C}}_{2,0}(i,j)+\frac{1}{2}({\tilde{C}}_{2,2}(i,j)+{\tilde{C}}_{2,-2}(i,j))\right)\\ & C_{y,y}(i,j)=-\frac{1}{\sqrt{3}}\left({\tilde{C}}_{0,0}(i,j)-\frac{1}{\sqrt{6}}{\tilde{C}}_{2,0}(i,j)-\frac{1}{2}({\tilde{C}}_{2,2}(i,j)+{\tilde{C}}_{2,-2}(i,j))\right)\\ & C_{z,z}(i,j)=-\frac{1}{\sqrt{3}}\left({\tilde{C}}_{0,0}(i,j)+\sqrt{\frac{2}{3}}{\tilde{C}}_{2,0}(i,j)\right).\end{array}$$

Full text

Introduction—Driven by the rapidly advancing research in artificial intelligence, many-body physics has embraced machine learning (ML) as a powerful tool to tackle non-trivial problems Carleo et al. (2019a). Applications include the use of neural networks for phase classification Carrasquilla and Melko (2017); van Nieuwenburg et al. (2017); Morningstar and Melko (2018); Beach et al. (2018); Ch’ng et al. (2017); Venderley et al. (2018); Casert et al. (2019); Liu and van Nieuwenburg (2018); Zhang et al. (2018), accelerating Monte Carlo algorithms Wang (2017); Liu et al. (2017a, b); Huang and Wang (2017); Bojesen (2018), and ML-based generative modeling of distributions in many-body physics Torlai and Melko (2016); Liu et al. (2017c); Mills and Tamblyn (2017); Carleo and Troyer (2017). The connection between the renormalization group and (deep) learning has also been highlighted Koch-Janusz and Ringel (2018); Li and Wang (2018); Mehta and Schwab (2014); Iso et al. (2018).

Quantum mechanical spin systems play a key role in the field of many-body physics. In Ref. Carleo and Troyer (2017), a particular class of artificial neural networks (ANN), namely restricted Boltzmann machines (RBM), is introduced as a variational ansatz for wave functions of many-body spin systems. The versatility of the ANN ansatz has been illustrated in studies of bosonic systems Saito (2017); Saito and Kato (2017), (chiral) topological states Glasser et al. (2018); Deng et al. (2017a); Zheng et al. (2018), frustrated systems Liang et al. (2018); Choo et al. (2019) and open systems Hartmann and Carleo (2019); Yoshioka and Hamazaki (2019); Nagy and Savona (2019); Vicentini et al. (2019). The RBM ansatz has been studied from the perspective of entanglement Deng et al. (2017b), and was shown to embody volume-law entanglement. In this light, the connections and differences between RBMs and matrix product states have been laid out Chen et al. (2018); Pastori et al. (2019).

The invariance of wave functions under symmetries, in particular $\textrm{SU}(2)$ symmetry, is important for applications such as quantum chemistry Szabó and Ostlund (1996); Gunst et al. (2019) and the description of spin liquids Zhou et al. (2017); Savary and Balents (2017). In quantum chemistry, wave functions are eigenfunctions of the total angular momentum operators $\widehat{J}^{2}$ and $\widehat{J}_{z}$, and spin liquid states do not break any symmetry of the Hamiltonian. The capacity of RBMs to capture long-range correlations and their independence of the problem’s geometry make them a prime candidate for these applications.

In this Letter, we introduce a methodology for constructing an RBM variational wave function that transforms as an irreducible representation of $\textrm{SU}(2)$. This wave function is designed to have a well-defined total angular momentum. Hence the method provides direct access to the construction of excited states. The challenge of imposing physical symmetries on the ANN ansatz has been addressed for finite abelian symmetry groups Choo et al. (2018), but this approach is not directly applicable to nonabelian symmetries. The proposed ansatz is not restricted to $\textrm{SU}(2)$ symmetry and can be applied to other nonabelian symmetries, as well as anyonic spin chains.

RBM wave functions—We use RBMs, which are energy-based generative neural networks, as a variational ansatz for quantum many-body spin systems. In the context of this work, the RBM models a distribution $\Psi(x_{1},...,x_{N})\equiv\Psi(\mathbf{x})$ of the variables $x_{i\in\{1,...,N\}}$, characterized by the energy function

[TABLE]

Here, $h_{j\in\{1,...,M\}}\in\{-1,1\}$ is a set of binary latent variables. The set $\mathcal{W}=\{w_{ij},a_{i},b_{j}\}$ are variational parameters: $w_{ij}$ are the weights connecting variables $x_{i}$ and $h_{j}$, and $a_{i}$ ($b_{j}$) are the biases of the physical (latent) variables $x_{i}$ ($h_{j}$). The ratio of the number of latent variables $M$ and the number of physical variables $N$, defined by $\alpha\equiv M/N$, is a measure of the complexity of the model.

The RBM was introduced in Ref. Carleo and Troyer (2017) as a variational ansatz for quantum many-body wave functions by modeling the probability amplitudes $\Psi(\mathbf{x})$ of the wave function $\ket{\Psi}=\sum_{\mathbf{x}}\Psi(\mathbf{x})\ket{\mathbf{x}}$ as the marginalized Boltzmann distribution

[TABLE]

This model is an extension of the RBM that is used to represent classical probability distributions as a marginalization over hidden variables of a Boltzmann distribution. In this work, $x_{i}$ is not necessarily a binary variable but remains discrete. We therefore replace $w_{ij}x_{i}$ with a set of general discrete functions $w_{ij}[x_{i}]$. This has the effect of transforming the variables $x_{i}$ non-linearly before constructing the RBM energy function of Eq. (1).

To find eigenstates of a Hamiltonian $\widehat{H}$, we use the variational principle to minimize the energy functional

[TABLE]

with respect to the parameters $\mathcal{W}$, for which we use the stochastic reconfiguration method Sorella (1998).

$\textrm{SU}(2)$* symmetry*—Hamiltonians with spin-rotation symmetry are omnipresent throughout the various domains of quantum many-body physics. Examples of $\textrm{SU}(2)$-symmetric quantum-mechanical spin chains are the antiferromagnetic Heisenberg model (AFH), the bilinear-biquadratic model, and the Majumdar-Ghosh model. In order to illustrate the potential of the proposed methodology, we focus on the AFH Hamiltonian

[TABLE]

with spin-$\nicefrac{{1}}{{2}}$ and spin-1 degrees of freedom in one dimension, which we study with open boundary conditions.

The main contribution of this work is the construction of a method to include $\textrm{SU}(2)$ symmetry in ANN wave functions. This procedure yields wave functions $\ket{J\,M_{J}}$ with well-defined total angular momentum $J$ and projection $M_{J}$. These wave functions belong to the subspace spanned by states with quantum numbers $J$ and $M_{J}$ of the full Hilbert space of the system. Part of the energy spectrum of the system can be uncovered by finding the variational minima in the subspaces with fixed quantum numbers $\ket{J\,M_{J}}$. For example, the gap of the AFH can be calculated after constructing the ground state $\ket{J=0\,M_{J}=0}$ and the first excited state $\ket{J=1\,M_{J}\in\{-1,0,1\}}$.

To construct wave functions $\ket{J\,M_{J}}$, we use Clebsch-Gordan coefficients to couple two spins, $\hat{\mathbf{s}}_{1}$ and $\hat{\mathbf{s}}_{2}$, to a single angular momentum degree of freedom $\hat{\mathbf{j}}_{2}=\hat{\mathbf{s}}_{1}+\hat{\mathbf{s}}_{2}$

[TABLE]

Given a system of $N$ spins $\hat{\mathbf{s}}_{i\in\{1,...,N\}}$, we construct states with total angular momentum $J$ by starting at the left side of the chain and using Eq. (5) to couple the first two spins to the angular momentum $j_{2}$. Next, $j_{2}$ is coupled to $s_{3}$. This process is repeated till reaching the end of the chain, resulting in a total angular momentum $J$. This is depicted in Fig. 1(a) with $y_{i}\equiv s_{i}$, $x_{i}\equiv j_{i}$ and $x_{0}=0$, $x_{N}=J$. Given the large amount of intermediate couplings, it is non-trivial and numerically challenging to transform the wave function in the coupled basis back into the $\hat{\mathbf{s}}_{i\in\{1,...,N\}}$ basis. As will become clear, however, observables of the studied quantum systems can be reliably and efficiently computed in the basis of coupled angular momenta. In the ansatz proposed in Ref. Carleo and Troyer (2017), the spin projections $m_{s_{i}}$ of $\hat{\mathbf{s}}_{i}$ are used as input values for the RBM in Eq. (2). Rather, we use the intermediate degrees of freedom $j_{k}$ as input (Fig. 1(b)), which produces the wave function

[TABLE]

where $\sum_{j_{k}}$ denotes a summation over all physically allowed configurations $j_{k\in\{1,...,N-1\}}$. Eq. (6) transforms as an irreducible representation of $\textrm{SU}(2)$, labeled by total angular momentum $J$, with dimension $2J+1$. For the states with $J=0$ (of which the ground state of the spin-$\nicefrac{{1}}{{2}}$ AFH is an example), the state is manifestly invariant under $\textrm{SU}(2)$ transformations, as the irreducible representation has dimension $1$. More information on this basis transformation can be found in the Supplemental Material (SM) SM . The above procedure can be readily extended to other nonabelian symmetries by decomposing the degrees of freedom in irreducible representations of the symmetry group and finding the equivalent of the Clebsch-Gordan coefficients in Eq. (5) to relate the irreducible representation of a system to those of its subsystems.

Spin-$\nicefrac{{1}}{{2}}$ AFH—Finding the ground state of the spin-$\nicefrac{{1}}{{2}}$ AFH in the coupled spin basis amounts to variationally minimizing the energy functional of Eq. (3) in the subspace defined by total angular momentum $J=0$. The accuracy of the variational wave function can be assessed by comparison to energies obtained with exact diagonalization (ED) $\Delta E_{0}=|(E_{0}-E_{0,exact})/E_{0,exact}|$, the magnitude of the variance of the Hamiltonian $\textrm{Var}(\widehat{H})=\langle\widehat{H}^{2}\rangle-\langle\widehat{H}\rangle^{2}$, and the weight $\varepsilon$ of excited states in the wave function. The latter can be found by writing $\ket{\Psi}=\sqrt{1-\varepsilon^{2}}\ket{\Psi_{GS}}+\varepsilon\ket{\Psi_{\perp}}$, where $\ket{\Psi_{GS}}$ is the exact ground state wave function and $\ket{\Psi_{\perp}}$ is a normalized superposition of states perpendicular to $\ket{\Psi_{GS}}$. The parameter $\varepsilon$ is a measure for the accuracy of the variational ground state as it measures the spurious content in $\ket{\Psi}$. An upper bound on $\varepsilon$ is given by the relation $\textrm{Var}(\widehat{H})\geq\varepsilon^{2}G^{2}$, with $G\equiv E_{1}-E_{0}$ the difference between the energies of the first excited state and the ground state.

In Fig. 2(a-c), we compare these convergence criteria for the ground states obtained with the RBM ansatz in the coupled basis and in the $s_{z}$-basis, as a function of the ratio $\alpha$. Both the relative energy error and $\textrm{Var}(\widehat{H})$ are systematically lower when using the coupled basis compared to the $s_{z}$-basis. For small systems, where ED is feasible, the parameter $\varepsilon$ can be determined exactly by computing the overlap $\braket{\Psi_{GS}}{\Psi}=\sqrt{1-\varepsilon}$. We obtain consistently lower values of $\varepsilon$ in the coupled basis. Related to this, we see that the spin-spin correlators are consistently better described, and inherently unbiased, in the coupled basis. This is described in the SM SM .

Also with the eye on gaining profound insight in the structure of the wave function, our methodology offers opportunities by studying the weight of the expansion coefficients in the ansatz of Eq. (6). We find that the basis state with the largest modulus has all pairs of neighboring spins coupled to a singlet. Next in importance are states with two neighboring triplets coupled to a singlet, on a background of singlets. More information on the structure of the wave function can be found in the SM SM .

The introduction of the coupled basis allows us to find the variational minimum of the energy functional of Eq. (3) in a subspace with specific $\ket{J\,M_{J}}$, which enables us to construct excited states. We demonstrate this by calculating the energy difference between the lowest lying eigenstate in the subspaces defined by $\ket{J=1\,M_{J}=0}$ and $\ket{J=0\,M_{J}=0}$. As the AFH is critical, the gap $G$ vanishes as $G\propto N^{-1}$ for $N\rightarrow\infty$. The gap as a function of system size is depicted in Fig. 3(a) and matches results obtained with ED or density matrix renormalization group (DMRG). The relative energy errors $\Delta E_{0}$ on the ground state range from $\mathcal{O}(10^{-5})$ for the smallest system sizes to $\mathcal{O}(10^{-4})$ for larger systems. The errors on the excited state energies $\Delta E_{1}$ are generally slightly larger.

Spin-1 AFH—Physically, the spin-1 AFH is inherently different from the spin-$\nicefrac{{1}}{{2}}$ AFH. Whereas the spin-$\nicefrac{{1}}{{2}}$ AFH is gapless in the thermodynamic limit, the spin-1 AFH has a fourfold degenerate ground state (consisting of a spin singlet and a spin triplet), above which a gap exists. The degeneracy of the ground state arises from the presence of effective spin-$\nicefrac{{1}}{{2}}$ degrees of freedom at the edges of the system. The interaction between these effective degrees of freedom is exponentially suppressed with system size, resulting in two free spin-$\nicefrac{{1}}{{2}}$ degrees of freedom in the thermodynamic limit White and Huse (1993). For finite systems, the ground state is non-degenerate and a spin singlet, in accordance with Marshall’s theorem Marshall (1955). The physical differences between the spin-$\nicefrac{{1}}{{2}}$ and spin-1 AFH make it interesting to investigate the representational ability of RBMs in both cases. Figs. 2(d-f) shows $\Delta E$, $\textrm{Var}(\widehat{H})$ and $\varepsilon$ as a function of the ratio $\alpha$ in the $s_{z}$-basis and the coupled basis. Across the whole range of $\alpha\in[0.5,4]$, the level of accuracy quantified by these measures is improved by at least an order of magnitude in the coupled basis as compared to the $s_{z}$-basis. These results are indicative for the effectiveness of the coupled basis. The structure of the wave function is similar to that of the spin-$\nicefrac{{1}}{{2}}$ Heisenberg model, and is described in detail in the SM SM .

Figure 3(b) shows the energy gap $G$ between the ground state $\ket{J=0\,M_{J}=0}$ and the first excited state $\ket{J=1\,M_{J}=0}$. For the energy gap, excellent agreement is reached between the RBM and ED/DMRG methods for different system sizes. The relative energy error on the ground state is $\mathcal{O}(10^{-6})$ for the smallest system sizes, and settles to $\mathcal{O}(10^{-5})$ for larger system sizes, while that of the excited state is $\mathcal{O}(10^{-5})$ to $\mathcal{O}(10^{-4})$. In Fig. 3(c), the energy gap of the spin-1 AFH with physical spin-$\nicefrac{{1}}{{2}}$ degrees of freedom on the edges is shown. The introduction of spin-$\nicefrac{{1}}{{2}}$ edges lifts the degeneracy of the ground state, and introduces a gap in the system, corresponding to the Haldane gap in the thermodynamic limit. Fig. 3(c) shows that the RBM ansatz in the coupled basis can represent gapped systems accurately.

Anyonic golden chain—As a novel application of the RBM ansatz, we turn to the simulation of anyonic systems. Anyons are defined as degrees of freedom that do not obey the mutual statistics of fermions or bosons Trebst et al. (2008). They arise in the fractional quantum Hall effect Arovas et al. (1984) and play a role in quantum computation Kitaev (2003); Nayak et al. (2008). Here, we study Fibonacci anyons, defined by one anyon type $\tau$ along with the trivial vacuum state $\mathbb{1}$. Anyonic many-particle systems are most easily described by fusing the different anyons. For Fibonacci anyons, the fusion rules are given by $\mathbb{1}\otimes\tau=\tau\otimes\mathbb{1}=\tau$ and $\tau\otimes\tau=\mathbb{1}\oplus\tau$. For a one-dimensional system, fusing the $N$ anyons can be done in a linear fashion from left to right (Fig. 1(a)), where $y_{i\in\{1,...,N\}}\equiv\tau$ and $x_{i\in\{1,...,N-1\}}\equiv\mathbb{1},\tau$. This construction is called a fusion tree. We choose the boundary conditions $x_{0}=x_{N}=\mathbb{1}$. Interacting anyons can be described on the level of the fused anyons. An example is the golden chain, which is reminiscent of the AFH model, defined by the Hamiltonian

[TABLE]

where $\widehat{\mathcal{P}}^{\mathbb{1}}_{i,i+1}$ is the projector on the vacuum fusion channel of the anyons with indices $i$ and $i+1$. The model defined by $\widehat{H}_{GC}$ is critical.

We exploit RBMs as a variational ansatz to find the ground state of $\widehat{H}_{GC}$. The results of the intermediate fusings are used as input for the RBM (Fig. 1(b)). The relative error on the ground-state energy $\Delta E_{0}$, and the variance of the Hamiltonian $\textrm{Var}(\widehat{H}_{GC})$ are shown in Fig. 4 for different values of $\alpha=M/N$. The relative energy error is below $\mathcal{O}(10^{-3})$ for all values of $\alpha$, reaching $\mathcal{O}(10^{-5})$ for $\alpha=4$. The variance of the Hamiltonian follows the same trend, ranging from $\mathcal{O}(10^{-3})$ to $\mathcal{O}(10^{-4})$.

Conclusion—In this paper, we have extended the variational class of artificial neural network states for spin systems to include nonabelian symmetries, in particular the $\textrm{SU}(2)$ spin-rotation symmetry. Hereto, we have formulated a restricted Boltzmann machine (RBM) ansatz using configurations labeled by intermediate spin coupling quantum numbers. Our numerical findings provide convincing support for using these quantities, which do not depend on the choice of local basis for the spins, as input variables to the RBM. Indeed, our ansatz obtains a higher precision in the ground-state energy for the same amount of variational freedom, compared to previous studies using the basis-dependent local spin projections as input variables. In doing so, the states we construct have well-defined total angular momentum $J$ and projection $M_{J}$ quantum numbers and can also be used to target the lowest-lying excited states of the system. Another application which makes use of the specific structure exploited by our ansatz are anyonic chains, for which our ansatz also accurately captures the ground state.

The approach presented here can be used to model spin liquid states, which are invariant under symmetries of the Hamiltonian. Likewise, our method could be adopted for the determination of wave functions in a quantum chemistry context. Both these applications are particularly suitable to be studied with RBMs, due to their potential to model long-range entanglement and correlations, and due to the geometrically independent way of modelling correlations in RBMs. Furthermore, our approach can be generalized for systems in more than one dimension, by defining an appropriate coupling scheme. These include coupling in the form of a graph, such as a tensor network Singh and Vidal (2012), or coupling in a tree-like fashion. We stress that our approach is independent of the form of the variational state and is applicable to other variational wave functions.

Appendix A Basis transformation

To find eigenstates of the $\textrm{SU}(2)$ symmetry group with restricted Boltzmann machines (RBM), the spins in the system are coupled to well-defined angular momentum eigenstates $\ket{J\,M_{J}}$ of the total system. Two spins $\hat{\mathbf{s}}_{1}$ and $\hat{\mathbf{s}}_{2}$ can be coupled to a total angular momentum $\hat{\mathbf{j}}_{2}=\hat{\mathbf{s}}_{1}+\hat{\mathbf{s}}_{2}$ via $3j$-symbols, or equivalently Clebsch-Gordan coefficients, as

[TABLE]

Eq. (8) relates product states of two basis states $\ket{s_{1}\,m_{s_{1}}}$ and $\ket{s_{2}\,m_{s_{2}}}$ of the Hilbert spaces $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ to states $\ket{s_{1}\,s_{2};j_{2}\,m_{j_{2}}}$ representing the tensor product space $\mathcal{H}_{1}\otimes\mathcal{H}_{2}$. The states $\ket{s_{1}\,s_{2};j_{2}\,m_{j_{2}}}$ are angular-momentum eigenstates, implying that $\hat{\mathbf{j}}_{2}^{2}\ket{s_{1}\,s_{2};j_{2}\,m_{j_{2}}}=j_{2}(j_{2}+1)\ket{s_{1}\,s_{2};j_{2}\,m_{j_{2}}}$ and $\hat{j}_{z_{2}}\ket{s_{1}\,s_{2};j_{2}\,m_{j_{2}}}=m_{j_{2}}\ket{s_{1}\,s_{2};j_{2}\,m_{j_{2}}}$. The orthonormality of the states $\ket{s_{1}\,s_{2};j_{2}\,m_{j_{2}}}$ in the tensor product space $\mathcal{H}_{1}\otimes\mathcal{H}_{2}$ follows readily from its definition and the use of the completeness relations of the tensor product basis

[TABLE]

The completeness of the coupled basis follows from the completeness of the product states and the fact that the coupled basis has the same number of orthonormal elements.

A system comprising $N$ spins can be coupled to a total angular momentum $\widehat{\mathbf{J}}=\hat{\mathbf{s}}_{1}+\hat{\mathbf{s}}_{2}+...+\hat{\mathbf{s}}_{N}$ by sequentially using the coupling rule of Eq. (8) until all the spins are coupled to a total angular momentum state $\ket{J\,M_{J}}$. The sequential coupling of spins can be performed using different schemes. We adopt a scheme that consists of coupling the spins $\hat{\mathbf{s}}_{i}$ in a linear fashion from left to right. Specifically, we start with an ancillary angular momentum state $\ket{j_{0}\,m_{j_{0}}}=\ket{0\,0}$ to which we couple $\ket{s_{1}\,m_{s_{1}}}$, resulting in $\ket{j_{1}=s_{1}\,m_{j_{1}}=m_{s_{1}}}$. Next, $\ket{s_{2}\,m_{s_{2}}}$ is coupled to $\ket{j_{1}\,m_{j_{1}}}$ resulting in $j_{2}\in\{0,1\}$. This is repeated until the end of the chain is reached, where we couple $\ket{s_{N}\,m_{s_{N}}}$ to $\ket{j_{N-1}\,m_{j_{N-1}}}$, resulting in $\ket{j_{N}\equiv J\,m_{j_{N}}\equiv M_{J}}$. This coupling scheme yields the set of intermediate total angular momenta $\{j_{1},j_{2},...,j_{N-1}\}$ as degrees of freedom. According to the angular momentum addition rules, the possible configurations of these intermediate total angular momenta have to fulfill the triangle inequalities $|j_{i-1}-s_{i}|\leq j_{i}\leq|j_{i-1}+s_{i}|$.

With this coupling scheme, the full basis transformation can be written as

[TABLE]

The expression of Eq. (10) is obtained by repeatedly using Eq. (8) to couple the intermediate angular momenta. The procedure is outlined in the above-mentioned coupling scheme and is graphically depicted in Fig. 1(a) of the main text.

From now on, we denote $\ket{j_{0}\,...\,j_{N-1};J\,M_{J}}\equiv\ket{s_{1}\,...\,s_{N}\,j_{0}\,...\,j_{N-1};J\,M_{J}}$, as the spin degrees of freedom $s_{1},...,s_{N}$ are fixed. The basis $\ket{j_{0}\,...\,j_{N-1};J\,M_{J}}$ forms an orthonormal basis, as can be seen from

[TABLE]

where we have used the orthogonality relation of $3j$-symbols

[TABLE]

from left to right. The completeness of the basis follows from the fact that, given the product states of two complete bases, the Clebsch-Gordan coefficients relate these states to a complete basis of the coupled system. Using the Clebsch-Gordan coefficients consecutively thus retains the completeness when all subsystems are described by a complete basis.

Appendix B F-move

To find the matrix elements of the antiferromagnetic Heisenberg Hamiltonian

[TABLE]

the basis is recoupled in such a way that the physical spins $\hat{\mathbf{s}}_{i}$ and $\hat{\mathbf{s}}_{i+1}$ that are subjected to the interaction $\hat{\mathbf{s}}_{i}\cdot\hat{\mathbf{s}}_{i+1}$ are coupled to a total angular momentum $\widehat{\mathbf{S}}=\hat{\mathbf{s}}_{i}+\hat{\mathbf{s}}_{i+1}$. This can be done by using an F-move, which is a recoupling using $6j$-symbols.

Graphically and mathematically, the F-move adopts the form of Fig. 5, where

[TABLE]

The term between curly brackets is a $6j$-symbol. This yields the basis transformation

[TABLE]

With $\hat{\mathbf{s}}_{i}\cdot\hat{\mathbf{s}}_{i+1}=\frac{1}{2}(S^{2}-s_{i}^{2}-s_{i+1}^{2})$ and the recoupling scheme of Eq. (15) one finds

[TABLE]

Using the recoupling defined in Eq. (15), the matrix elements are found in the basis labeled by states $\ket{j_{0}\,j_{1}\,...\,j_{N-1};J\,M_{J}}$

[TABLE]

As can be seen in Eq. (17), the matrix elements of the operator $\hat{\mathbf{s}}_{i}\cdot\hat{\mathbf{s}}_{i+1}$ are independent of the values of $j_{k\notin\{i-1,i,i+1\}}$, are diagonal in $j_{k\in\{i-1,i+1\}}$ and are not diagonal in $j_{i}$.

Appendix C $\textrm{SU}(2)$-symmetry

The $\textrm{SU}(2)$ symmetry group is a continuous symmetry group. For spin degrees of freedom, it is represented by the unitary operators

[TABLE]

where $\hat{s}_{x,y,z}$ are the generators of $\textrm{SU}(2)$ and $\bm{\theta}$ is a normalised vector in three-dimensional space. For spin-$\nicefrac{{1}}{{2}}$, $\hat{s}_{x,y,z}=\hat{\sigma}_{x,y,z}/2$, where $\hat{\sigma}_{x,y,z}$ are the two-dimensional Pauli-matrices

[TABLE]

For spin-$1$, $\hat{s}_{x,y,z}$ are the three-dimensional matrices

[TABLE]

The representation of the Hilbert space in terms of basis states with a total angular momentum forms an irreducible representation with respect to the $\textrm{SU}(2)$ symmetry, leaving the total angular momentum invariant and transforming the angular momentum projection degrees of freedom. For a recoupling of spins using Clebsch-Gordan coefficients

[TABLE]

acting with the unitary operator $\widehat{U}(\bm{\theta})=e^{i\hat{\mathbf{j}}\cdot\bm{\theta}}=e^{i(\hat{\mathbf{s}}_{1}+\hat{\mathbf{s}}_{2})\cdot\bm{\theta}}$ on both sides of the equation yields

[TABLE]

This means that the coupling which transforms $\ket{s_{1}\,m_{s_{1}},s_{2}\,m_{s_{2}}}$ in $\ket{s_{1}\,s_{2};j\,m_{j}}$ also transforms $e^{i(\hat{\mathbf{s}}_{1}+\hat{\mathbf{s}}_{2})\cdot\bm{\theta}}\ket{s_{1}\,m_{s_{1}},s_{2}\,m_{s_{2}}}$ in $e^{i\hat{\mathbf{j}}\cdot\bm{\theta}}\ket{s_{1}\,s_{2};j\,m_{j}}$. Using this result sequentially on the coupling used in this work, we see that the coupling which transforms $\ket{s_{1}\,m_{s_{1}};...;s_{N}\,m_{s_{N}}}$ in $\ket{s_{1}\,...\,s_{N}\,j_{0}\,j_{1}\,...\,j_{N-1};J\,m_{J}}$ also transforms $e^{i(\hat{\mathbf{s}}_{1}+...+\hat{\mathbf{s}}_{N})\cdot\bm{\theta}}\ket{s_{1}\,m_{s_{1}};...;s_{N}\,m_{s_{N}}}$ in $e^{i\hat{\mathbf{J}}\cdot\bm{\theta}}\ket{s_{1}\,...\,s_{N}\,j_{1}\,...\,j_{N-1};J\,m_{J}}$. In the specific case where we couple to the angular momentum state $\ket{J=0,M_{J}=0}$, $e^{i\hat{\mathbf{0}}\cdot\bm{\theta}}=\widehat{\mathbb{1}}$, which does indeed leave the coupled state invariant. This proves that the variational ground-state wave function is manifestly invariant under $\textrm{SU}(2)$ transformations since every basis state is itself invariant under $\textrm{SU}(2)$ transformations.

Appendix D Basis cut-off

The basis transformation used in this paper introduces the intermediate angular momenta $j_{i\in\{1,...,N-1\}}$ as degrees of freedom. In the case of spin-$\nicefrac{{1}}{{2}}$, these can take on values $j_{i}\in\{0,\nicefrac{{1}}{{2}},1,\nicefrac{{3}}{{2}},2,...\}$, while in the case of spin-$1$, they can take on values $j_{i}\in\{0,1,2,3,...\}$. As explained earlier, the intermediate degrees of freedom need to satisfy the triangle relations $\left|j_{i-1}-s_{i}\right|\leq j_{i}\leq\left|j_{i-1}+s_{i}\right|$. These constraints induce a maximal value $j_{max}$ for the degrees of freedom, which is proportional with system size $N$. Because the RBM parameters explicitly depend on the input values, the number of weights also increases linearly with $j_{max}$.

The eigenstate with minimal eigenvalue of the individual terms of the Hamiltonian is the state where the two spins involved couple to a singlet. Hence, one expects physically that the antiferromagnetic Heisenberg Hamiltonian favors the coupling of spins to a small total angular momentum. This motivates the introduction of a cut-off $j_{cut}$ for the maximal value of the intermediate angular momenta $j_{i}$. In Fig. 6, the effect of introducing the cut-off $j_{cut}$ is inspected by comparing the relative energy error of the ground state energy of the AFH with and without this constraint, for different system sizes. The results are obtained with exact diagonalization. We see that the relative energy error decreases rapidly and settles to machine precision at $j_{cut}=3$ for the spin-$\nicefrac{{1}}{{2}}$ model and $j_{cut}=4$ for the spin-1 model. Moreover, the dependence on system size is negligible, even for the largest systems, indicating that the cut-off is inherent to the physics of the system.

Appendix E Wave function structure

A central question in variational many-body state optimization is the level of accuracy of observables other than the energy (which is used as the optimization criterion). More information can be obtained by looking at how close the individual expansion coefficients of the chosen basis functions lie to the exact expansion coefficients. Hence, we compare the expansion coefficients of the RBM model to those of exact diagonalization in the coupled basis. An interesting follow-up question is which configurations of intermediate spins have the highest importance in the basis expansion of the ground state, where we define the importance as the modulus of the expansion coefficient. We denote $\Psi_{j}$ as the expansion coefficient with the $j$-th largest modulus. Because the RBM wave function is not normalized, we compare the square modulus of the expansion coefficients relative to that with the largest square modulus, i.e. $|\Psi_{j}|^{2}/|\Psi_{0}|^{2}$. Fig. 7 shows the relative importance of the expansion coefficients $|\Psi_{j}|^{2}/|\Psi_{0}|^{2}$ as a function of $j$ for both the variational wave function and the wave function obtained with exact diagonalization. For $|\Psi_{j}|^{2}/|\Psi_{0}|^{2}\gtrsim 10^{-6}$ the expansion coefficients coincide, while for $|\Psi_{j}|^{2}/|\Psi_{0}|^{2}\lesssim 10^{-6}$ they start to diverge slightly. The left column of Fig. 8 shows from top to bottom the $11$ most important configurations of the spin-$\nicefrac{{1}}{{2}}$ AFH in the coupled basis. Using the recoupling of spins, we can provide an interpretation to the depicted configurations. The most important configuration can be recoupled to the case where, from left to right, every two neighbouring spins are coupled to a singlet. Note that this is a particular case of the resonating valence bond state. The next $10$ configurations can all be recoupled to a background consisting of singlets (as in the first configuration), but where two neighbouring singlets are excited to two neighbouring triplets, which couple together to a singlet. The results for the spin-1 AFH are shown in the right column of Fig. 8, and are qualitatively similar to the spin-$\nicefrac{{1}}{{2}}$ case.

Appendix F Correlations

To assess the advantage of imposing symmetries in RBM wave functions, we compute the spin-spin correlation functions for the wave functions discussed in the main text. The spin-spin correlation functions are defined as

[TABLE]

with $\beta,\gamma\in\{x,y,z\}$ the different directions of the spin operator $\hat{\mathbf{s}}$. For the Heisenberg model, the components of the correlation function with $\beta=\gamma$ should be equal due to the inherent $\textrm{SU}(2)$-symmetry. To prove that, in the coupled basis, this principle holds by construction, we make use of the Wigner-Eckart theorem

[TABLE]

The Wigner-Eckart theorem states that the overlap of a $\tilde{T}_{l,m_{l}}$ spherical tensor operator of rank $l$ and component $m_{l}$ between two angular momentum eigenstates $\ket{J\,m_{J}}$ and $\ket{J^{\prime}\,m^{\prime}_{J}}$ is equal to the Clebsch-Gordan coefficient $\braket{J^{\prime}\,m^{\prime}_{J};l\,m_{l}}{J\,m_{J}}$ multiplied with a reduced matrix element $\braket{J^{\prime}}{\widehat{T}_{l}}{J}$, independent of the angular momentum projections. As one can see from Eq. (24), using states with $J=J^{\prime}=0$, which is the case for the ground state in the coupled basis, all expectation values of spherical operators with $l\neq 0$ are identically zero.

We write the $9$ Cartesian components of the correlation functions in terms of the operators in the spherical basis via a linear transformation. Here, we use the definitions of Ref. Emsly (1985). In particular, the spherical operator with $l=0$ and $m_{l}=0$ is defined as

[TABLE]

The Cartesian operators with $\beta=\gamma$ are related to the spherical operators by

[TABLE]

Using the Wigner-Eckart theorem for the Cartesian operators of Eq. (26), we find that, as only the spherical component with $(l=0,m_{l}=0)$ contribute, all three Cartesian components in Eq. (26) are equal to a third of the scalar product $\mathbf{\hat{s}}_{i}\cdot\mathbf{\hat{s}_{j}}$.

In Fig. 9, we show the correlation functions and the absolute error of the correlation functions with respect to the exact values for both the spin-$\nicefrac{{1}}{{2}}$ (Fig. 9 (a-b)) and spin-$1$ (Fig. 9 (c-d)) Heisenberg model, as described in the main text. In particular, we show the correlation functions $C(i,j)$ as a function of distance $|i-j|$ for the lowest ($\alpha=0.5$) and for the highest ($\alpha=4$) number of hidden units. For the $s_{z}$-basis, we show the longitudinal ($C^{s_{z}}_{\parallel}(i,j)=C^{s_{z}}_{z,z}(i,j)$) and transverse ($C^{s_{z}}_{\perp}(i,j)=(C^{s_{z}}_{x,x}(i,j)+C^{s_{z}}_{y,y}(i,j))/2$) correlation functions separately. For the coupled basis, we show $C^{\textrm{SU(2)}}(i,j)=C^{\textrm{SU(2)}}_{x,x}(i,j)=C^{\textrm{SU(2)}}_{y,y}(i,j)=C^{\textrm{SU(2)}}_{z,z}(i,j)$. For the correlations in the spin-$\nicefrac{{1}}{{2}}$ Heisenberg model, one can notice that the RBM with $\alpha=0.5$ in the $s_{z}$-basis breaks the symmetry between the longitudinal and transverse components and is biased for the correlations in the longitudinal direction. This is a direct consequence of using the $s_{z}$-basis as input states for the variational wavefunction. For the RBM in the coupled basis, the correlations lie close to the exact values, even for $\alpha=0.5$. For $\alpha=4$, both RBMs perform well, with absolute errors of $\mathcal{O}(10^{-5})$. However, the inequality between the transverse and longitudinal components in the $s_{z}$-basis persists. For the correlations with $\alpha=0.5$ in the spin-$1$ Heisenberg model, the deviations are much larger in the $s_{z}$-basis, also for the longitudinal component. While the scalar product of spins are reasonably well described, the individual components are largely over- or underestimated. By construction, this is not the case in the coupled basis, which is further evidence for the power of using coupled irreducible representations of the symmetry groups of the model as basis states for the variational model.

In all our experiments, the dependence of the error on the correlation functions on distance is weak. This is in contrast with the correlation functions in the transverse field Ising model, which were recently studied in ref.Collura et al. (2019). However, as the models and system sizes are different, it is difficult to compare both results.

Appendix G Hyperparameters

Several hyperparameters have to be chosen for the variational optimization of RBMs. These hyperparameters are the learning rate $\eta$ of the optimization algorithm (stochastic gradient descent), the diagonal shift $d$ of the stochastic reconfiguration algorithm, and the standard deviation $\sigma$ of the normal distribution (centered at 0) from which the parameters of the RBM are initialized. In order to find the set of hyperparameters which enables the energy to reach its lowest values, we draw a number (order 100) of random combinations of hyperparameters, and use these to optimize RBM wave functions. Typically, we sample $\eta$ uniformly in the range $\eta\in[0.01,0.2]$, $\log(d)$ uniformly in the range $\log(d)\in[-3,-1]$, and $\sigma$ uniformly in the range $\sigma\in[0.01,0.1]$. We found that the energy after optimization is largely independent of $d$ and $\sigma$, but depends strongly on the learning rate $\eta$. For the experiments in the coupled basis, we reached the lowest energies with $\eta\approx 0.1$, while for those in the $s_{z}$-basis, we used learning rates $\eta\in[0.01,0.1]$. Finally, we used parallel tempering Swendsen and Wang (1986) for the optimization of the spin-1 AFH. The parallel tempering algorithm uses a set of Monte Carlo chains (replicas) that all sample the probability distribution raised to a power between [math] and $1$. As a result, the chains with lower power sample a smoother distribution than those with higher power. By mixing configurations between neighbouring chains (according to a Monte Carlo step), the chain with power $1$ (the physical distribution which is used to measure observables) samples a richer region of the phase space. The number of replicas in the parallel tempering algorithm was 50 and were chosen with uniformly separated powers.