Continuous-time Monte Carlo Renormalization Group

In recent years, continuous-time quantum Monte Carlo (QMC) has become the standard to simulate sign-free lattice quantum spin systems Rieger and Kawashima (1999); Blöte and Deng (2002); Mishchenko et al. (2000); Greitemann and Pollet (2018). In continuous-time QMC, one adopts a path-integral representation, in which the time dimension is represented by a genuine continuous line of length equal to the inverse temperature $\beta$, whereas the space dimension is represented by a discrete lattice. For systems that exhibit Lorentz invariance at criticality, the path integral represents a statistical field theory, in which isotropy between the time and space directions emerges at large distance. This occurs in systems whose low-lying elementary energy excitations depend linearly on the magnitude of the momentum of the excitation:

Here, $E_{0}$ is the ground state energy, and ${\bf k}$ is the momentum of a low-lying energy eigenstate of energy $E({\bf k})$. The constant of proportionality, $v_{s}$, is called the sound velocity and is a non-universal constant specific to a system. The sound velocity is the scale factor connecting time and length scales of a system. When Eq. 1 happens, the large-distance behavior of a lattice spin system is described by a massless quantum field theory, invariant under dilational coordinate transformations. The invariance under dilational transformation underlies the success of renormalization group (RG) theory Wilson (1971) in treating critical phenomena, as done in the real space Monte Carlo Renormalization Group (MCRG) approach Swendsen (1979), which is often used to estimate critical couplings and exponents of spin systems. When implemented on a discrete lattice, the dilational transformation is limited to integral scale factors. However, non-integral scale factors are desirable in many cases. For example, anisotropic coarse graining can lead to a fixed point Hamiltonian isotropic in the space and time directions. Then, the ratio of space and time lengths in the anisotropic coarse graining is the sound velocity. The value of the sound velocity is generally a real number, requiring that the coarse-graining either along the space or the time direction be continuous. We show that this can be realized, within the MCRG framework, for quantum spin models on a lattice by adopting a continuous imaginary time path integral representation. In this way, Lorentz invariance is demonstrated numerically and accurate estimates of the sound velocity are obtained. In addition, continuous time MCRG allows us to compute the lattice version of the energy stress tensor of the underlying field theory.

Conformal invariance is very powerful in two-dimensions (2D), due to the infinite dimensionality of the local conformal algebra Belavin et al. (1984). Conformal field theory (CFT) yields finite sizing scaling predictions of physical observables, such as the energy Affleck (1986); Blöte et al. (1986) and the entanglement entropy Calabrese and Cardy (2004). Sound velocity and the energy-stress tensor are often parameters in these predictions. Thus, assuming the validity of these predictions, the sound velocity and the energy-stress tensor may be obtained by fitting observables in a numerical simulation against the CFT predictions. This has been carried out in many studies, e.g. von Gehlen et al. (1986); Ma and He (2019). With continuous-time MCRG, we determine these quantities from their defining expressions, without recourse to CFT results. This allows a much easier generalization to three dimensions (3D), where CFT predictions are less available. In this paper, we illustrate the idea of continuous-time MCRG mostly with examples in 2D, showing that already in 2D our approach leads to estimates of the sound velocity that are more accurate than those obtained with alternative methods, such as directly computing the energy spectrum or fitting numerical observables against CFT predictions. We also provide a 3D example, by reporting a calculation of the sound velocity for the quantum Ising model in two space dimensions, a system for which, to the best of our knowledge, the sound velocity cannot be obtained with other means. A more detailed study of systems in 3D is deferred to future works.

The paper is organized as follows. In Sec. 2, we use a diagrammatic expansion to obtain the path integral representation of the partition function of the $Q$-state Potts model, the system that we use here as an example to illustrate the methodology. In Sec. III, we coarse grain the time direction and explain the MCRG procedure. In Sec. IV, we use continuous-time MCRG to compute the sound velocity. In Sec. V, we interpret the coarse-graining along the time direction as a continuous coordinate transformation and discuss its connection with the energy-stress tensor. In Sec. VI, we report our conclusions.

For concreteness, we illustrate the method with the $Q$-state Potts model having Hamiltonian Sólyom and Pfeuty (1981):

where the system is in a $d$-dimensional hypercubic lattice with length $L$. $i$ and $j$ label different lattice sites, and $\langle i,j\rangle$ denotes a nearest neighbor bond. The operators $\hat{\Omega}_{i}$ and $\hat{M}_{i}$ act on the $Q$ states of the local Hilbert space at site $i$, which we label by $|0\rangle_{i},...,|s\rangle_{i},...|Q-1\rangle_{i}$. In this local basis, the $\hat{\Omega}_{i}$ is a diagonal matrix such that $\hat{\Omega}_{i}|s\rangle_{i}=\omega^{s}|s\rangle_{i}$, where $\omega=e^{i2\pi/Q}$ and $s=0,\cdots,Q-1$. $\hat{M}_{i}$ performs a cyclic permutation: $|0\rangle_{i}\rightarrow|Q-1\rangle_{i},|1\rangle_{i}\rightarrow|0\rangle_{i},\cdots,|Q-1\rangle_{i}\rightarrow|Q-2\rangle_{i}$, and acts as a transverse-field. When $d=1$, this model is self-dual and a quantum phase transition occurs at the critical coupling $g_{c}=1$ for all $Q$ Sólyom and Pfeuty (1981). When $Q\leq 4$, the transition is continuous and is described by a CFT Baxter (1973); von Gehlen et al. (1986). When $Q>4$, the transition is first-order and has a finite latent heat at $g_{c}$ Baxter (1973).

To derive a path-integral representation of the system partition function $Z=\text{Tr}(e^{-\beta\hat{H}_{Q}})$, one takes

and considers the partition function in the interaction picture

where $\hat{T}$ is the time-ordering operator in imaginary time $\tau$, and $\hat{H}_{1}(\tau)=e^{\tau\hat{H}_{0}}\hat{H}_{1}e^{-\tau\hat{H}_{0}}$ is the $\hat{H}_{1}$, in the interaction picture. Eq. 4 can be written as a diagrammatic expansion in the following way:

Here the states $|\bm{\sigma}\rangle=\otimes_{i}|\sigma\rangle_{i}$ form a basis in the Hilbert space. Each index $i_{1},\cdots,i_{n}$ runs from $1$ to $L^{d}$. $\hat{H}_{0}$ is diagonal in the $|\bm{\sigma}\rangle$ basis: $\hat{H}_{0}|\bm{\sigma}\rangle=Q\sum_{\langle i,j\rangle}\delta_{\sigma_{i},\sigma_{j}}|\bm{\sigma}\rangle\equiv h_{0}(\bm{\sigma})|\bm{\sigma}\rangle$. Eq. 5 suggests the following Monte Carlo (MC) scheme to sample the partition function. For each $i=1,\cdots,L^{d}$ and each $\tau\in[0,\beta]$, a Potts spin, $\sigma_{i}(\tau)$, ranging from $0$ to $Q-1$, is assigned to an MC configuration. As the state $|\bm{\sigma}\rangle$ is propagated in imaginary time, spin flips can happen at any lattice site and at any time. Let $\tau_{l}$ and $i_{l}$ denote the $l$th flip time and its associated lattice site. Here $l$ could be equal to $1,2,\cdots,$ or $n$. In addition, let $\tau_{l}^{-}$ and $\tau_{l}^{+}$ respectively denote the time immediately before and after the flip time $\tau_{l}$. If a spin flip occurs at $\tau_{l}$ on site $i_{l}$, $\sigma_{i_{l}}(\tau_{l}^{-})$ will be made to switch to any $\sigma_{i_{l}}(\tau_{l}^{+})$ different from $\sigma_{i_{l}}(\tau_{l}^{-})$ by the action of $\sum_{k=1}^{Q-1}\hat{M}_{i_{l}}^{k}$. In Eq. 5, the earliest spin flip occurs at $\tau_{1}$ on site $i_{1}$; the second one occurs at $\tau_{2}$ on site $i_{2}$, etc.. The total number of spin flips, $n$, contributes a weight $g^{n}$ to the sampling of the diagrammatic expansion. In addition, the weight includes factors equal to $e^{-(\tau_{l+1}-\tau_{l})h_{0}(\bm{\sigma}(\tau_{l}^{+}))}$ between two consecutive spin flips at $\tau_{l}$ ad $\tau_{l+1}$. Finally, the periodicity of the trace requires $\bm{\sigma}(\beta)=\bm{\sigma}(0)$, which in turn implies that $n$ should be even. Thus, MC sampling does not have a sign problem even if $g$ is negative.

The partition function in Eq. 5 is given by a sum of terms (diagrams) that entail summation over discrete variables and integration over continuous ones. The contribution of the different terms, which are associated to the weights detailed above, is evaluated stochastically with a MC algorithm that follows the protocols discussed in Refs. Rieger and Kawashima (1999); Blöte and Deng (2002); Mishchenko et al. (2000). For the $Q$-state Potts model diagrammatic MC can use a continuous time cluster algorithm Blöte and Deng (2002), based on the Wolff algorithm Wolff (1989), which significantly reduces equilibration time. We will use both local and cluster MC algorithms in the following.

Eq. 5 indicates that the thermodynamics of a $d$-dimensional quantum Potts model is described by a statistical field theory in $(d+1)$ dimensions, where on each lattice site $i$ there is a worldline of length $\beta$ described by the function $\sigma_{i}(\tau)$. We can coarse grain this worldline by placing $P$ lattice points along the time direction through the majority-rule. That is, we partition the worldline into $P$ intervals: $[0,\frac{P}{\beta}],[\frac{P}{\beta},\frac{2\beta}{P}],\cdots,[\frac{(P-1)\beta}{P},\beta]$. In each MC configuration, to each interval, we assign the Potts spin which appears most often on that interval. By discretizing time in this way we represent each worldline with $P$ discrete Potts spins and we end up with a ($d$+1)-dimensional hypercubic lattice that hosts $PL^{d}$ Potts spin. Each MC snapshot corresponds to a configuration of those spins. The probability distribution of the spin configurations on the discrete lattice is not known explicitly, but can be sampled by coarse-graining the configurations generated in the diagrammatic MC simulation Blöte and Deng (2002). We can then perform MCRG on the $(d+1)$-dimensional hypercubic lattice. We designate the $n$-th RG iteration with a subscript $(n)$, where $n=0,1,2,..$. In the 0-th iteration the spin configuration, $\bm{\sigma}^{(0)}$, is the one obtained from coarse graining the time direction of diagrammatic MC. The probability distribution of $\bm{\sigma}^{(0)}$ is described by the (unknown) lattice Hamiltonian $H^{(0)}$. The subsequent levels of coarse graining are generated by successive block spin RG transformations and are characterized by spin configurations $\bm{\sigma}^{(n)}$ and Hamiltonians $H^{(n)}$.

In all of the above coarse-graining transformations we use short-ranged coupling terms $S_{\alpha}(\bm{\sigma})$ to parametrize the probability distribution $P(\bm{\sigma})$ of the spin configurations $\bm{\sigma}$:

with coupling constants $K_{\alpha}$. The terms $S_{\alpha}(\bm{\sigma})$ include nearest neighbor (nn), next nearest neighbor (nnn), smallest plaquette ($\square$), etc., interactions.

On the hypercubic lattice, a dilation transformation with scale factor $b$ is realized through a block spin transformation by a conditional probability, $T(\bm{\mu}|\bm{\sigma})$, of the renormalized spin $\bm{\mu}$ given unrenormalized spins $\bm{\sigma}$. For example, in the majority-rule coarse-graining that we use, $\bm{\sigma}$ is partitioned into hypercubic blocks with side $b$. When there is a unique spin $\sigma_{\text{max}}$ that appears the most often in the block, the renormalized spin $\mu$ for that block is assigned to $\sigma_{\text{max}}$ with probability one. When $p$ $\sigma$-spins tie for the most appearances in the block, $\mu$ is assigned to one of these $p$ spins with probability $\frac{1}{p}$. Assuming the lattice is large enough, the spin configurations may be iteratively coarse-grained $n$ times, from $\bm{\sigma}^{(0)}$ to $\bm{\sigma}^{(n)}$, giving rise to a renormalized Hamiltionian at the $n$th RG level:

where $T^{(n)}(\bm{\sigma}^{(n)}|\bm{\sigma}^{(0)})$ is the conditional probability of $\bm{\sigma}^{(n)}$ given $\bm{\sigma}^{(0)}$. $T^{(n)}(\bm{\sigma}^{(n)}|\bm{\sigma}^{(0)})$ implements coarse-graining from $\bm{\sigma}^{(0)}$ to $\bm{\sigma}^{(n)}$, realizing a dilation transformation with scale factor $b^{n}$. It is obtained by iterating the single-level coarse-graining $n$ times:

In principle, an infinite number of coupling terms $S_{\alpha}(\bm{\sigma}^{(n)})$ is required to exactly parametrize $H^{(n)}$. In practice, this is not possible, and we only use a few coupling terms, such as those listed in Table 1. Given a set of coupling terms, we use Variational Monte Carlo Renormalization Group (VMCRG) Wu and Car (2017) to compute the corresponding coupling constants $K^{(n)}_{\alpha}$. VMCRG allows us to compute the renormalized coupling constants in a way that greatly alleviates the critical slowing down and has very little finite-size effect Wu and Car (2019). Its implementation details are explained in Wu and Car (2017, 2019); Wu (2019). Here by finite-size effect, we mean the effect of different system size, $L$, with the same RG level, $n$. Our result, however, necessarily depends on what $n$ one uses to approach the fixed-point Hamiltonian. Using only a finite number of coupling terms in the Hamiltonian in Eq. 6 introduces a truncation errors. However, because the truncation scheme respects isotropy, i.e. the truncation of an isotropic Hamiltonian is still isotropic, we can estimate the sound velocity without truncation errors.

The sound velocity $v_{s}$ is an important property of a quantum system whose low lying excited states show linear momentum dispersion. In particular, this quantity is required to compare the predictions from CFT with observables of a lattice model. When $d=1$, one can compute the sound velocity with finite size scaling of the ground state energy or entanglement entropy, assuming validity of the CFT prediction Affleck (1986); Blöte et al. (1986); Calabrese and Cardy (2004), or one can directly compute the excitation spectrum of the system. The latter calculation can be done by exact diagonalization of the system Hamiltonian, which is limited to small lattice sizes, or it can be done with density matrix renormalization group (DMRG) techniques Schollwöck (2011), which introduce truncation errors due to finite bond dimension. When $d>1$, neither method works reliably, and one has to resort to QMC. In fact, the sound velocity has been calculated with continuous-time QMC by looking for a scale factor $v_{s}$ such that the coorelation function $C(x,{\bf v}_{s}\tau)$ becomes isotropic along the time and space directions at large distances Prokof’ev et al. (1998); Deng and Blöte (2002). In the following, we will use directly RG to compute the sound velocity and show that VMCRG can deal with rather large lattice sizes, up to at least $L=256$, leading to very accurate estimates of the sound velocity.

To compute the sound velocity, we perform a dilation transformation with $b=2$ using the majority rule. In the VMCRG calculation, we take couplings along space and time directions to be independent, as the system is necessarily anisotropic between space and time. However, in the presence of Lorentz invariance, isotropy between space and time is recovered at large distances up to a scale factor $v_{s}$. One can thus vary $\frac{P}{\beta}$ until the couplings along the space direction and those along the time direction become equal for a large $n$. The $\frac{P}{\beta}$ so determined is the sound velocity, $v_{s}$.

As one changes the parameter $\frac{P}{\beta}$ in the zeroth RG iteration, one also changes the fixed-point Hamiltonian reached by the RG procedure, as shown, for example, in Table 1. Since dilational transformations are isotropic, there is a line of fixed-point Hamiltonians reflecting the different extent of anisotropy in the system Cardy (1996). A change of $\frac{P}{\beta}$ generates a movement along this line of fixed-point Hamiltonians. In fact, fixing $P$, the change $\beta\rightarrow\beta+\delta\beta$ induces a coordinate transformation: $x_{0}\rightarrow x^{\prime}_{0}=(1-\frac{\delta\beta}{\beta})x_{0},x_{1}\rightarrow x^{\prime}_{1}=x_{1}$, where $x_{0}$ and $x_{1}$ are time and space coordinates, respectively. Here we have taken the coordinate transformation to be passive, i.e. ${\bf x}=(x_{0},x_{1})$ and ${\bf x}^{\prime}=(x^{\prime}_{0},x^{\prime}_{1})$ denote the number of lattice spacings needed to describe the same physical length, before and after the transformation. Thus, a time dilation generates a change in the system Hamiltonian. In field theory, the response of the system Hamiltonian to a generic coordinate transformation, $x^{\mu}\rightarrow x^{\prime\mu}=x^{\mu}+\epsilon^{\mu}({\bf x})$, is described by the energy-stress tensor, $T^{\mu\nu}$, defined by

where $D$ is the space-time dimension of the system. As $\beta$ is conjugate to $\hat{H}$ in the action, we identify $T_{00}$ as the energy operator in the path integral.

To appreciate the novelty brought by VMCRG in this context, let us consider, for example, the two-dimensional classical Ising model with the Hamiltonian

where $\langle i,j\rangle_{0}$ and $\langle i,j\rangle_{1}$ are nearest neighbor spins along the $x_{0}$ and the $x_{1}$ direction, respectively. The system is isotropic and critical when $K_{0}=K_{1}=K_{c}=\text{arcsinh}(1)/2=0.4407\cdots$ Schultz et al. (1964). An infinitesimal change in the coupling constant, $K_{0}=K_{c}-\delta J$ and $K_{1}=K_{c}+\delta J$, turns on anisotropy yet the system still maintains its criticality. That is, the deviation from the isotropic Hamiltonian,

generates a length scale transformation $x_{0}\rightarrow x^{\prime}_{0}=(1-\delta\lambda)x_{0}$ and $x_{1}\rightarrow x^{\prime}_{1}=(1+\delta\lambda)x_{1}$, where $\delta\lambda=\frac{1}{\gamma}\delta J$ with an unknown proportionality constant $\gamma$. Continuous-time VMCRG provides a way to determine $\gamma$ directly, and in that sense, it is a ruler of anisotropy.

To determine $\gamma$, we invoke the universality of the fixed-point Hamiltonians. Let $H^{*}(\bm{\mu})$ be the fixed-point Hamiltonian that VMCRG eventually reaches, starting from the critical $d=1$ TFIM with $\frac{P}{\beta}=v_{s}$ and from the critical isotropic $d=2$ classical Ising model. In practice, we approximate $H^{*}(\bm{\mu})$ with $H^{(n)}(\bm{\mu})$ for some large $n$. For the TFIM, the change in the action $\beta\hat{H}\rightarrow(\beta+\delta\beta)\hat{H}$ generates a change in the fixed-point Hamiltonian $\delta H^{(n)}(\bm{\mu})=-\sum_{\alpha}\frac{\partial K^{(n)}_{\alpha}}{\partial\beta}S_{\alpha}(\bm{\mu})\delta\beta$ with an anisotropy of the extent $\delta(\frac{x_{1}}{x_{0}})=\frac{x^{\prime}_{1}}{x^{\prime}_{0}}-\frac{x_{1}}{x_{0}}=\frac{\delta\beta}{\beta}\frac{x_{1}}{x_{0}}$. For the classical $d=2$ Ising model, the change in the unrenormalized Hamiltonian $H_{\text{Ising}}\rightarrow H_{\text{Ising}}+\sum_{m,n}A_{m,n}\delta J$ generates a change in the fixed-point Hamiltonian $\delta H^{(n)}(\bm{\mu})=-\sum_{\alpha}\frac{\partial K^{(n)}_{\alpha}}{\partial J}S_{\alpha}(\bm{\mu})\delta J$, with an anisotropy of the extent $\delta(\frac{x_{1}}{x_{0}})=2\delta\lambda\frac{x_{1}}{x_{0}}$. The $\delta H^{(n)}(\bm{\mu})$ for the TFIM and for the $d=2$ classical Ising model should be multiples of each other, because the line of fixed-point Hamiltonians is universal. In particular, when they are equal, the anisotropies that they represent should coincide. This means that

for all $\alpha$. The Jacobians of the RG transformation, $\frac{\partial K_{\alpha}^{(n)}}{\partial J}$ and $\frac{\partial K_{\alpha}^{(n)}}{\partial\beta}$, can be readily computed with VMCRG Wu and Car (2017), where the first Jacobian is calculated for the TFIM, and the second one for the classical Ising model. For the operator $A(\bm{\sigma})$ defined in Eq. 13, $\gamma$ is analytically known and is $\frac{\sqrt{2}}{2}=0.7071\cdots$. A VMCRG calculation using Eq. 14 with $n=4$ gives $\gamma=0.708\pm 0.001$, so the ruler works.

With the coordinate transformation $x_{0}\rightarrow x^{\prime}_{0}=(1-\delta\lambda)x_{0}$ and $x_{1}\rightarrow x^{\prime}_{1}=(1+\delta\lambda)x_{1}$, a part of the energy-stress tensor can now be read off from Eq. 11:

We take the lattice spacing to be 1, and $\sum_{m,n}$ is equivalent to $\int d^{2}x$. This gives

where, in 2D, $T$ and $\bar{T}$ are respectively the holomorphic and the antiholomorphic component of the energy-stress tensor, and are defined as $T=\frac{1}{4}(T_{00}-2iT_{01}-T_{11})$ and $\bar{T}=\frac{1}{4}(T_{00}+2iT_{01}-T_{11})$. While the argument is developed for the Ising model, it also generalizes to other 2D systems.

A consequence of Eq. 16 is that one obtains a prediction of the finite-size dependence of $\langle A\rangle$, due to CFT. For example, if one simulates a critical system infinitely long along the $x_{0}$ direction but periodic of size $L$ along $x_{1}$, CFT predicts that $\langle T\rangle=\langle\bar{T}\rangle=-(\frac{2\pi}{L})^{2}\frac{c}{24}$, and thus $\langle A\rangle=-\frac{1}{\gamma\pi}(\frac{2\pi}{L})^{2}\frac{c}{12}$, where $c$ is the central charge of the underlying CFT. This prediction on $A$ has been verified in Bastiaansen and Knops (1998).

In this paper, we have shown how to perform MCRG with continuous-time Monte Carlo simulations, and demonstrated that space-time isotropy is explicitly recovered at large distances. This yields a practical method to determine the sound velocity and the energy-stress tensor from their defining expressions. This should allow generalizations to systems in three dimensions, which could be studied in the future.