Quantum Simulation/Quantum Monte-Carlo

Resource type: this resource contains a lecture or lecture notes.

Quantum-classical mapping

We have seen in the previous chapter that exact diagonalization is an impossible task for more than a few particles. In quantum Monte-Carlo simulations, the goal is to avoid considering the full Hilbert space, but randomly sample the most relevant degrees of freedom and try to extract the quantities of interest such as the magnetization by averaging over a stochastic process. To pursue this goal, we first need to establish a framework in which we can interpet quantum-mechanical observables in terms of (classical) probabilities. This process is called a "quantum-classical mapping" and allows us to reformulate quantum many-body problems in terms of models from classical statistical mechanic, albeit in higher dimensions.

Suppose we wish to calculate the partition function $Z=\mathrm {Tr} \{\exp(-\beta H)\}$ of the transverse Ising model, as knowledge of the partition function allows us to calculate all thermodynamic quanties that may be of interest. Specifically, we have

${\begin{aligned}Z&=\mathrm {Tr} \left\{\exp \left[-\beta \left(g\sum \limits _{i}\sigma _{x}^{(i)}-\sum \limits _{i}\sigma _{z}^{(i)}\sigma _{z}^{(i+1)}\right)\right]\right\}=\mathrm {Tr} \left\{\exp \left(-{\frac {\beta g}{N_{y}}}\sum \limits _{i}\sigma _{x}^{(i)}+{\frac {\beta }{N_{y}}}\sum \limits _{i}\sigma _{z}^{(i)}\sigma _{z}^{(i+1)}\right)^{N_{y}}\right\}\\&=\lim \limits _{N_{y}\to \infty }\mathrm {Tr} \left\{\left[\exp \left(-{\frac {\beta g}{N_{y}}}\sum \limits _{i}\sigma _{x}^{(i)}\right)\exp \left({\frac {\beta }{N_{y}}}\sum \limits _{i}\sigma _{z}^{(i)}\sigma _{z}^{(i+1)}\right)\right]^{N_{y}}\right\},\end{aligned}}$

where in the last line we have used the Suzuki-Trotter expansion

$\exp \left[{\frac {1}{N}}(A+B)\right]=\exp \left({\frac {A}{N}}\right)\exp \left({\frac {B}{N}}\right)+O(1/N^{2}),$

which can be proved using the Baker-Campbell-Hausdorff formula. Using the same expansion, we can replace the exponential of the product by a product of the exponentials, i.e.,

$Z=\lim \limits _{N_{y}\to \infty }\mathrm {Tr} \left\{\exp \left(-{\frac {\beta g}{N_{y}}}\sum \limits _{i}\sigma _{x}^{(i)}\right)^{N_{y}}\exp \left({\frac {\beta }{N_{y}}}\sum \limits _{i}\sigma _{z}^{(i)}\sigma _{z}^{(i+1)}\right)^{N_{y}}\right\}.$

(1)

The exponentiation to the power of $N_{y}$ can be written as a product, and we can insert $N_{y}-1$ identity operators according to

$A^{N_{y}}=\prod \limits _{i=1}^{N_{y}}A=A|i_{1}\rangle \langle i_{1}|A|i_{2}\rangle \langle i_{2}|A\cdots A|i_{N_{y}-1}\rangle \langle i_{N_{y}-1}|A.$

Let us now look more closely at the product involving the spin-flip operators $\sigma _{x}$, where we will encounter terms of the form

$\langle i_{j}|\exp \left(a\sigma _{x}^{(i)}\right)|i_{j+1}\rangle =\langle i_{j}|\left[\cosh(a)+\sigma _{x}^{(i)}\sinh(a)\right]|i_{j+1}\rangle$.

The crucial part of the quantum-classical mapping is to interpret the partition function of the one-dimensional chain containing $N$ spins as the partition function of a corresponding two-dimensional spin model containing $N\times N_{y}$ spin . In this interpretation, we can rewrite the spin-flip operators in terms of an Ising interaction in the $y$ direction (plus a constant energy shift),

$\langle i_{j}|\left[\cosh(a)+\sigma _{x}^{(i)}\sinh(a)\right]|i_{j+1}\rangle ={\frac {1}{2}}[\exp(-a)\langle i_{j}|\sigma _{z}^{(i,j)}\sigma _{z}^{(i,j+1)}|i_{j+1}\rangle +\exp(a)].$

We now want to cast these terms back into an exponential form,

${\frac {1}{2}}[\exp(-a)\langle i_{j}|\sigma _{z}^{(i,j)}\sigma _{z}^{(i,j+1)}|i_{j+1}\rangle +\exp(a)]=\Lambda \exp \left(\gamma \sigma _{z}^{(i,j)}\sigma _{z}^{(i,j+1)}\right)$

where we find for the coefficents $\Lambda$ and $\gamma$

${\begin{aligned}\Lambda &={\sqrt {\sinh(a)\cosh(a)}}\\\gamma &=-{\frac {1}{2}}\log \tanh(a).\end{aligned}}$

We can now insert this expression back into the partition function Eq.(1) and carry out the $N_{y}$ multiplications. The boundary conditions for the Ising interaction along the $N_{y}$ direction are fixed by the final trace operation; as the trace is implemented by multiplying $\langle i_{N_{y}}|$ from the left and $|i_{N_{y}}\rangle$ from the right, we have periodic boundary conditions. In total, we obtain

$Z=\Lambda ^{NN_{y}}\mathrm {Tr} \left\{\exp \left(\gamma \sum \limits _{i=1}^{N}\sum \limits _{j=1}^{N_{y}}\sigma _{z}^{(i,j)}\sigma _{z}^{(i,j+1)}+{\frac {\beta }{N_{y}}}\sum \limits _{i=1}^{N}\sum \limits _{j=1}^{N_{y}}\sigma _{z}^{(i,j)}\sigma _{z}^{(i+1,j)}\right)\right\}.$

The constant prefactor in the partition function is irrelevant as it will drop out when calculating thermodynamic observables. Consequently, we can identify a corresponding two-dimensional classical Ising model with anisotropic interactions which reproduces the same thermodynamics as the one-dimensional quantum Ising model in a transverse field. The Hamiltonian for the classical model is given by

$H_{cl}=-{\frac {N_{y}\gamma }{\beta }}\sum \limits _{i=1}^{N}\sum \limits _{j=1}^{N_{y}}\sigma _{z}^{(i,j)}\sigma _{z}^{(i,j+1)}-\sum \limits _{i=1}^{N}\sum \limits _{j=1}^{N_{y}}\sigma _{z}^{(i,j)}\sigma _{z}^{(i+1,j)}.$

Note, however, that the classical temperature $\beta _{cl}=\beta /N_{y}$ is different from the quantum temperature $\beta$. Nevertheless, we can now proceed to calculate the thermodynamic properties of the quantum model by performing classical Monte-Carlo simulations.

Metropolis algorithm

When trying to evaluate thermodynamic observables for a classical spin model, we still find ourselves in considerable difficulties as also the classical configuration space grows exponentially with system size. However, we are not really interested in a solution that incorporates all microscopic details, but rather we want to obtain information about macroscopic observables. So, we will be fine with any microscopic description of the model of interest, as long as it gets the macroscopic statistics right. Here, the goal is to find a microscopic description which can be efficiently (i.e., using resources that only grow polynomially with system size) implemented on a computer.

The most famous method for the Monte-Carlo simulation of statistical mechanics models is the Metropolis algorithm . Let us first state the basic steps of the algorithm for the Ising model and then analyze it in more detail.

1. Pick an arbitrary initial state (e.g., all spins polarized) and compute its energy $E$.
2. Flip a random spin and calculate the energy of the new configuration $E'$
3. If $E'<E$, always accept the new configuration.
4. If $E'>E$, accept the new configuration with probability $\exp(-\beta [E'-E])$.
5. Continue at step 2 until the macroscopic observables (averaged over a fixed number of steps) are equilibrated.

Two-dimensional classical Ising model at $\beta =1$.

Two-dimensional classical Ising model at $\beta =\beta _{c}=\log(1+{\sqrt {2}})/2$.

To evaluate the algorithm, let us consider two configurations $x$ and $x'$ with energies $E$ and $E'$,respectively. The probalities for these configurations are denoted by $p(x)$, and $p(x')$. Let us assume that $E<E'$, so the transition probability satisfy $p(x'\to x)=1$ and $p(x'\to x)=\exp(-\beta [E'-E])$ for the inverse process. In equilibrium, the system will satisfy detailed balance (absence of currents), i.e,

$p(x)p(x\to x')=p(x')p(x'\to x),$

or equivalently

${\frac {p(x')}{p(x)}}={\frac {p(x\to x')}{p(x'\to x)}}=\exp(-\beta [E'-E]),$

which reproduces the Boltzmann distribution and thus gives rise to the correct thermodynamic behavior. The convergence of the Metropolis algorithm can be improved by also including unphysical processes that flip a large domain of spins at once (so-called "cluster updates") , as constructing these processes from individual spin flips will consume a lot of time. Fig. 1 shows results of Monte-Carlo simulations for two different values of $\beta$ for the isotropic two-dimensional Ising model.

From classical to quantum phase transitions

Let us now come back to the anisotropic Ising model that we obtained using the quantum-classical mapping. If we play around with $N_{y}$ and $\beta _{cl}$, we find that the phase transition between the paramagnet and the ferromagnet occurs on a critical line that is given by

$\sinh(2\beta _{cl}^{c})\sinh \left(2\beta _{cl}^{c}{\frac {N_{y}\gamma }{\beta }}\right)=1,$

where $\beta _{cl}^{c}$ is the critical inverse temperature . Remarkably, the condition $N_{y}\to \infty$ at a fixed classical temperature $\beta _{c}=\beta /N_{y}$ means that the corresponding quantum phase transition takes place at $\beta =\infty$, i.e., at zero temperature. At finite temperature, the extension in the $N_{y}$ direction will be finite and the classical model belongs to the universality class of the one-dimensional Ising model, which does not exhibit a phase transition. If we substitute the definitions for the coupling constants into the equation for the critical line, we obtain

${\frac {\sinh(2\beta _{cl}^{c})}{\sinh(2g_{c}\beta _{cl}^{c})}},$

which immediately gives us the critical transverse field $g_{c}=1$, in agreement with the exact solution of the quantum model. From our knowledge of the classical 2D Ising model, we can also immediate extract the scaling of observables close to the critical point, e.g., for the magnetization

$m=m_{0}(g-g_{c})^{1/8},$

as the critical exponents for the 1D quantum model are identical to the ones of the classical 2D model.

The sign problem

While quantum Monte-Carlo methods are indeed very powerful for tackling a large variety of many-body problems, they are also constrained by inherent limitations. One possible problem we can run into is that the partition function of the classical model cannot be computed efficiently. For example, consider the classical Hamiltonian

$H=\sum \limits _{\langle i,j\rangle }J_{ij}\sigma _{z}^{(i)}\sigma _{z}^{(j)}$

on a three-dimensional cubic lattice, where the nearest-neighbor interactions $J_{ij}$ are randomly chosen from the values $\{-1,0,+1\}$. This model exhibits glassy behavior and finding its ground state is known to be NP-complete , i.e., it is widely believed that there is no efficient way to simulate this model. On the other hand, we can safely assume that precisely due to the glassy properties that makes the model hard to study, the ground state is irrelevant as a physical state as the system will take exponentially long to reach it even if we put it into heat bath at zero temperature. A more subtle issue can arise as the result of the quantum-classical mapping, which is known as the "sign problem".

To be explicit, let us assume that during our quantum-classical mapping procedure, we encounter the following term in the Hamiltonian:

$H=J\sigma _{+}^{(i)}\sigma _{-}^{(i+1)}+\mathrm {h.c.} ,$

where we have used the spin-flip operators $\sigma _{+}=|\uparrow \rangle \langle \downarrow |$ and $\sigma _{-}=\sigma _{+}^{\dagger }=|\downarrow \rangle \langle \uparrow |$. Within the quantum-classical mapping, we will encounter terms of the form

$\langle i_{j}|\exp \left(a\sigma _{+}^{(i)}\sigma _{-}^{(i+1)}\right)|i_{j+1}\rangle ={\frac {1}{2}}[\cosh(a)+1+(\cosh(a)-1)\sigma _{z}^{(i)}\sigma _{z}^{(i+1)}+\sinh(a)\sigma _{+}^{i}\sigma _{-}^{i+1}].$

The constant term and the Ising interaction are unproblematic and can be cast into the classical partition function as in the case of the transverse Ising model. Whether we can do the same for the spin-flip term (which will ultimately result in a four-body interaction for the classical Ising spins), depends on the sign of $J$. For the ferromagnetic case, $J<0$, the coefficient $a=-\beta J/N$ will be positive and we can turn this interaction into a term of the form $\Lambda \exp(-\beta _{cl}J_{cl}H_{cl})$ with $\Lambda$ being positive so the interpretation in terms of a classical partition function remains valid. On the other hand, for an antiferromagnetic interaction with $J>0$, we find $\Lambda$ to be negative and therefore we do not obtain the equivalent of a classical partition function. In some cases, we can work around this by using a unitary transformation that maps onto a model that does not exhibit this problem. For instance on a bipartite lattice, we can fix this problem by flipping spins on one of the sublattices by transforming $\sigma _{z}\to -\sigma _{z}$, turning the antiferromagnetic interaction into a ferromagnetic one. However, this does not work in general, nevertheless there is still a way to perform classical Monte-Carlo sampling. When sampling over our classical configuration space, the problematic contributions correspond to negative probabilities, i.e., when we compute the value of an observable $A$,

$\langle A\rangle ={\frac {\mathrm {Tr} \{A\exp(-\beta H)\}}{\mathrm {Tr} \{\exp(-\beta H)\}}}={\frac {\sum \limits _{i}A_{i}p_{i}}{\sum \limits _{i}p_{i}}},$

we find some of the $p_{i}$ to be negative. This can be resolved by taking the absolute value of these probabilities and dividing by the expectation value of the sign, i.e.,

$\langle A\rangle ={\frac {\sum \limits _{i}A_{i}\operatorname {sgn} p_{i}|p_{i}|/\sum _{i}|p_{i}|}{\sum \limits _{i}\operatorname {sgn} p_{i}|p_{i}|/\sum _{i}|p_{i}|}}\equiv {\frac {\langle As\rangle _{|p|}}{\langle s\rangle _{|p|}}}.$

However, this modification comes at a hefty price: as the uncertainty of the denominator will typically increase exponentially with the size of the system , we will no longer have an efficient algorithm.

Despite these difficulties, it may sometimes still be more favorable to perform Monte-Carlo simulation than performing exact diagonalization, but this certainly depends on details of the model at hand. Finally, the sign problem is not limited to models of the type outlined above. As a rough guidance, frustrated spin models and fermionic models tend to have a sign problem, while bosonic models tend to be free of it.