Download Quantum Simulations in Physics and more Lecture notes Physics in PDF only on Docsity!
Lecture Notes on Quantum Simulation
Hendrik Weimer∗ Institut f¨Germanyur Theoretische Physik, Leibniz Universit¨at Hannover, Appelstraße 2, 30167 Hannover,
These lecture notes were created for a graduate-level course on quantum simula- tion taught at Leibniz University Hannover in 2013. The first part of the course discusses various state of the art methods for the numerical description of many- body quantum systems. In particular, I explain successful applications and inherent limitations due to the exponential complexity of the many-body problem. In the second part of the course, I show how using highly controllable quantum system such as ultracold atoms will provide a way to overcome the limitations of classical simulation methods. I discuss several theoretical and experimental achievements and outline the road for future developments.
Contents
Chapter 1: Exact Diagonalization 3 1.1. The many-body problem 3 1.2. Exact and not-so-exact diagonalization 4 1.3. The power method 5 1.4. Lanczos algorithm 5
Chapter 2: Quantum Monte-Carlo 8 2.1. Quantum-classical mapping 8 2.2. Metropolis algorithm 10 ∗Electronic address: [email protected]
These lecture notes are available under the terms of the Creative Commons Attribution-ShareAlike 4. International license. Sharing of this work must contain attribution to the author and the arXiv identifier.
arXiv:1704.07260v1 [quant-ph] 24 Apr 2017
- 2.3. From classical to quantum phase transitions
- 2.4. The sign problem
- Chapter 3: Density-Matrix Renormalization Group
- 3.1. Variational method
- 3.2. Reduced density matrices
- 3.3. The DMRG algorithm
- 3.4. Matrix product states
- 3.5. Entanglement and the area law
- Chapter 4: Analog quantum simulation
- 4.1. Ultracold atoms
- 4.2. Optical lattices
- 4.3. The Bose-Hubbard model
- 4.4. Fermi-Hubbard model
- Chapter 5: Digital quantum simulation
- 5.1. Quantum logic gates
- 5.2. Digital simulation procedure
- 5.3. Implementation based on Rydberg atoms
- Acknowledgements
spin in a magnetic field. Consequently, the ground state is given by
|φg〉 =
∏^ N
i=
√^1
2 (| ↑〉i^ − | ↓〉i).^ (1.3)
In the limit of zero external field, g → 0, we just have to minimize the interaction energy, which leads to two distinct fully polarized ground states, |φ(1) g 〉 = | ↓↓↓ · · · 〉 and |φ(2) g 〉 = | ↑↑↑ · · · 〉. Note that the Hamiltonian does not distinguish between up and down spins (a so-called Z 2 symmetry), but the two possible ground states do. This symmetry breaking is a manifestation of the system being in two different quantum phases for g 1 (ferromagnet) and g 1 (paramagnet). From the exact solution, it is known that the quantum phase transition occurs at a critical transverse field of gc = 1.
1.2 Exact and not-so-exact diagonalization
The term “exact diagonalization” is often used in a slightly misleading manner. In general, to find the eigenvalues of a d-dimensional Hamiltonian, one has to find the roots to the characteristic polynomial of degree d, for which in general no exact solution can be found for d > 4. Of course, we can still hope to numerically approximate the eigenvalues to an arbitrary degree, but the fact that we have to work with computers operating with fixed precision numbers makes this endeavour substantially more complicated. Keeping that aside, finding the eigenvalues by solving the characteristic polynomial is a bad idea, as finding roots of high-degree polynomials is a numerically tricky task [2]. In fact, one of the most powerful methods for finding the roots of such a polynomial is to generate a matrix that has the same characteristic polynomial and find its eigenvalues using a different algorithm! A much better strategy is to find a unitary (or orthogonal, if all matrix elements of the Hamiltonian are real) transformation that makes the Hamiltonian diagonal, i.e.,
H → U †HU (1.4) The general strategy is to construct the matrix U in an iterative way,
H → U 1 † HU 1 → U 2 † U 1 † HU 1 U 2 → · · · (1.5) until the matrix becomes diagonal. The columns of U = U 1 U 2 U 3 · · · then contains the eigenvectors of H. There are many different algorithms for actually performing the diago-
nalization, and it is usually a good idea to resort to existing libraries (such as LAPACK) for this task. However, if we are only interested in finding the low-energy eigenvalues of a particular Hamiltonian, we can make a few simplifications that will allow for much faster computations.
1.3 The power method
Initially, we pick a random state |φ 0 〉 which has a very small but finite overlap with the ground state, i.e., 〈φ 0 |φg〉 6 = 0. Then, we repeatedly multiply the Hamiltonian with this initial state and normalize the result,
|φn+1〉 = N H|φn〉, (1.6)
where N is the normalization operation. This method will eventually converge to the eigen- vector with the largest absolute eigenvalue, so by subtracting a constant energy from the Hamiltonian, we can always ensure that this will be the ground state. The key advantage is that the matrix-vector multiplications occuring at each iteration can be implemented very fast: in most cases (as for the transverse Ising model), for each column the number of nonzero entries in such a sparse Hamiltonian are much smaller than the dimension of the Hilbert space (here: N vs. 2N^ ).
1.4 Lanczos algorithm
The power method can be readily improved by using not only a single state during each iteration, but employing a larger set of states which will be extended until convergence is reached. This procedure is known as Lanczos algorithm and is implemented as follows [3]:
- Pick a random state |φ 0 〉 as in the Power method
- Construct a second state |φ 1 〉 according to |φ 1 〉 = H|φ 0 〉 − 〈φ〈^0 φ| 0 H|φ|φ 0 〉^0 〉|φ 0 〉 (1.7)
- Starting with n = 2, recursively construct an orthogonal set of states given by
|φn+1〉 = H|φn〉 − an|φn〉 − b^2 n|φn− 1 〉, (1.8)
1.5 Some remarks on complexity
Let us briefly consider the difficulty inherent in exact diagonalization studies. On each lattice site, we have 2 degrees of freedom, refering to a spin pointing either up or down. However, for N spins, we have 2N^ possible spin configuration. The consequences of this exponential scaling cannot be underestimated: if we want to store the vector corresponding to the ground state in a computer using double precision, we would need 8 TB of memory for N = 40, and for N = 300, the number of basis states already exceeds the number of atoms in the universe! From these considerations, we see that exact diagonalization can only work in the limit of small system sizes.
References
[1] S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, 1999). [2] W. H. Press, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in C (Cambridge University Press, Cambridge, 1992). [3] E. Dagotto, Rev. Mod. Phys. 66 , 763 (1994).
CHAPTER 2: QUANTUM MONTE-CARLO
2.1 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 = Tr{exp(−β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
Z = Tr
exp
[
−β
g
i
σ x(i )−
i
σ z(i )σ z(i+1)
)]}
= Tr
exp
− (^) Nβgy
i
σ( xi )+ (^) Nβy
i
σ z(i )σ z(i+1)
)Ny^ (2.2)
= (^) Nlimy →∞ Tr
[
exp
− (^) Nβgy
i
σ( xi)
exp
β Ny
i
σ( zi )σ( zi+1)
)]Ny^ ,^ (2.3)
where in the last line we have used the Suzuki-Trotter expansion
exp
[ 1
N (A^ +^ B)
]
= exp
( A
N
exp
( B
N
+ O(1/N 2 ), (2.4)
which can be proved using the Baker-Campbell-Hausdorff formula. Using the same ex- pansion, we can replace the exponential of the product by a product of the exponentials, i.e.,
Z = (^) Nlimy →∞ Tr
exp
− βgNy
i
σ( xi)
)Ny exp
β Ny
i
σ( zi )σ( zi+1)
)Ny^ .^ (2.5)
The exponentiation to the power of Ny can be written as a product, and we can insert Ny − 1 identity operators according to
ANy^ =
∏^ Ny i=
A = A|i 1 〉〈i 1 |A|i 2 〉〈i 2 |A · · · A|iNy − 1 〉〈iNy − 1 |A. (2.6)
Note, however, that the classical temperature βcl = β/Ny is different from the quantum temperature β. Nevertheless, we can now proceed to calculate the thermodynamic properties of the quantum model by performing classical Monte-Carlo simulations.
2.2 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 ex- ponentially with system size. However, we are not really interested in a solution that incor- porates 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 [2]. Let us first state the basic steps of the algorithm for the Ising model and then analyze it in more detail.
- Pick an arbitrary initial state (e.g., all spins polarized) and compute its energy E.
- Flip a random spin and calculate the energy of the new configuration E′
- If E′^ <; E, always accept the new configuration.
- If E′^ > E, accept the new configuration with probability exp(−β[E′^ − E]).
- Continue at step 2 until the macroscopic observables (averaged over a fixed number of steps) are equilibrated. 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′^ → x) = 1 and p(x′^ → x) = exp(−β[E′^ − E]) for the inverse process. In equilibrium, the system will satisfy detailed balance (absence of currents), i.e,
p(x)p(x → x′) = p(x′)p(x′^ → x), (2.14)
FIG. 2.1 Two-dimensional classical Ising model at β = 1 (left) and at β = βc = log(1 + √2)/ 2 (right) (taken from [3]).
or equivalently p(x′) p(x) =^
p(x → x′) p(x′^ → x) = exp(−β[E
′ − E]), (2.15)
which reproduces the Boltzmann distribution and thus gives rise to the correct thermo- dynamic 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”) [4], 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 β for the isotropic two-dimensional Ising model.
on a three-dimensional cubic lattice, where the nearest-neighbor interactions Jij are ran- domly chosen from the values {− 1 , 0 , +1}. This model exhibits glassy behavior and finding its ground state is known to be NP-complete [5], 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 irrel- evant 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σ+(i) σ −(i+1) + h.c., (2.20)
where we have used the spin-flip operators σ+ = | ↑〉〈↓ | and σ− = σ† + = | ↓〉〈↑ |. Within the quantum-classical mapping, we will encounter terms of the form
〈ij | exp
aσ+(i) σ −(i+1)
|ij+1〉 =^12 [cosh(a)+1+(cosh(a)−1)σ( zi )σ( zi +1)+sinh(a)σ +iσ −i+1 ]. (2.21)
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 = −βJ/N will be positive and we can turn this interaction into a term of the form Λ exp(−βclJclHcl) with Λ 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 Λ 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 σz → −σ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,
〈A〉 = Tr Tr{{Aexp(^ exp(−−βHβH)})} =
i ∑^ Aipi i^ pi
we find some of the pi 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.,
〈A〉 =
i
Ai sgn pi|pi|/ ∑ i |pi| ∑ i^ sgn^ pi|pi|/^
i |pi|^
≡ 〈As 〈s〉〉|p| |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 [6], 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.
References
[1] G. G. Batrouni and R. T. Scalettar, in Ultracold Gases and Quantum Information, edited by C. Miniatura, L.-C. Kwek, M. Ducloy, B. Grmaud, B.-G. Englert, L. Cugliandolo, A. Ekert, and K. K. Phua (Oxford University Press, 2011). [2] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21 , 1087 (1953). [3] AlterVista, http://commons.wikimedia.org/wiki/File:Ising beta 1.gif and http://commons.wikimedia.org/wiki/File:Ising betaC.gif. [4] W. Krauth, in New optimization algorithms in physics, edited by A. K. Hartmann and H. Rieger (Wiley-VCH ; John Wiley, Weinheim; Chichester, 2004). [5] F. Barahona, J. Phys. A 15 , 3241 (1982). [6] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94 , 170201 (2005).
in the context of open quantum systems. Generally, we can think of a system decribed by a density matrix ρ as a statistical mixture of pure quantum states |ψi〉, i.e.,
ρ =
i
pi|ψi〉〈ψi|, (3.4)
where pi can be interpreted as the probability to find the system in the pure quantum state |ψi〉. Density matrices have unit trace and are positive-semidefinite, and their time evolution is governed by the von Neumann equation, d dtρ^ =^ −^
i ℏ[H, ρ].^ (3.5) Density matrices play an important role when looking at a smaller subsystem embedded in an environment. Consider a bipartite system composed of the subsystems A and B, with pure states being represented by
|ψ〉 =
ij
cij |i〉A|j〉B ≡
ij
cij |ij〉. (3.6)
Then, we can define the reduced density matrix of A as the partial trace over B, given by
ρA = TrB {ρ} =
i,i′
j
〈ij|ρ|i′j〉|i〉〈i′|. (3.7)
3.3 The DMRG algorithm
The key element of DMRG is to think of the many-body system of interest being composed of a system of size l, attached to an environment of the size L − l. Suppose we have already found a set of D states {|φS 〉}, which allows us to give a good approximation to both the Hamiltonian and its ground state for this particular system size. Then, we can increase the system size to l + 1 by the following procedure [3]:
- Form a new system S′^ of size l + 1 by combining the D states describing the system with the full Hilbert space describing the additional site. For spin 1/2, we have a new set of states according to {|φSl ↑〉, |φSl ↓〉}.
- Build a superblock of size 2l + 2 by combining the enlarged system with the enlarged environment. The Hilbert space of this dimension will be 4D^2 in the case of spins. If the Hamiltonian is reflection symmetric (i.e, left and right are identical), the basis states for the system and the environment will be the same.
- Find the ground state |ψ〉 of the Hamiltonian of this superblock by exact diagonaliza- tion.
- Compute the reduced density matrix for the enlarged system
ρ′ S = Tr′ E {|ψ〉〈ψ|}. (3.8)
- Diagonalize the reduced density matrix and keep only the eigenvectors corresponding to the D largest eigenvalues. This set of states {|ψSl+1〉} will serve as the input for the next iteration step.
- Express the Hamiltonian for system size l + 1 in this new basis.
- Continue this iteration until sufficient convergence of observables (e.g, of the ground state energy) is obtained. This procedure is commonly referred to as “infinite-system DMRG” [3]. In practice, one can often achieve significantly better results by accounting for finite size effects. This can be done in a relatively straightforward way by keeping track of the system and the environment separately. Once the infinite-system algorithm has reached the desired size, one can successively increase the system at the expense of the environment. When the environment reaches a minimum size, the procedure is reversed and the environment is increased at the expense of the system. Multiple sweeps of this kind can be performed until the ground state energy has converged to an even better value. During each step, we can also use the estimate for the ground state obtained during the previous step as the initial state of our exact diagonalization procedure, which will result in a drastically improved performance [3]. In general, the error of the DMRG procedure will be given by the truncation error ε, which is the weight of all the states that have been dropped during the DMRG step. Typical values are ε ∼ 10 −^10 for D ∼ 100. As in the case of exact diagonalization, DMRG is not restricted to ground states, but due to the constraints imposed by the convergence of the exact diagonalization step, it works best for low-energy eigenvalues. It is also possible to compute the time evolution of a many- body system using DMRG, for instance by combining it with the Runge-Kutta scheme for the numerical integration of ordinary differential equations [4]. This is done by using the
where M forms the largest block that can be solved exactly, i.e, dim(M ) = D. Then, we can interpret the coefficients ψmσl/ 2 σl/2+1nas the variational parameters according to which the energy is being minimized. Further information about the relationship between DMRG and matrix product states can be found in a more recent review article [6].
3.5 Entanglement and the area law
Having seen that DMRG implements a variational method based upon matrix product states, it is important to pose the question how generic these states will be in terms of ground states of many-body Hamiltonians. A key concept in answering this question is entanglement, i.e., nonclassical correlations between different parts of the system. Consider the two-spin state |ψ〉 = √^12 (| ↑↑〉 + | ↓↓〉, (3.18)
which can be represented by the density matrix
ρ =
^.^ (3.19)
As can be seen from the diagonal elements, this state possesses classical correlations, but its off-diagonal contributions are what is crucial here. It is instructive to look at the reduced density matrix of a single spin,
ρ 1 = Tr 2 {ρ} =
1 /^2
This state is proportional to the identity, so it is invariant under all unitary transformations, i.e., U ρU †^ = ρ. Consequently, despite the two-spin state being pure, its local density matrix of a single spin is completely random and does not carry any information at all! We can quantify this behavior using the von Neumann entropy, defined as
S = −Tr{ρ log ρ}, (3.21)
which vanishes for pure states and reaches its maximum log d, with d being the Hilbert space
dimension, for the “maximally mixed state”,
ρ =
1 /d 1 /d
...
which for a single spin is precisely the state we have obtained above. The von Neumann entropy is invariant under unitary transformations and therefore does not change under Hamiltonian dynamics. Additionally, it is subadditive with respect to its subsystems, i.e.,
S(ρAB ) ≤ S(ρA) + S(ρB ), (3.23)
where ρA and ρB are the density matrices of the individual subsystems A and B [7]. Conse- quently, by looking only at parts of the full system, we can only obtain partial information, as seen in the example above. Remarkably, for pure states of the full system, the entropies S(ρA) and S(ρB ) are identical. This can be seen by looking at the Schmidt decomposition [8] of the wave function of the full system,
|ψ〉 =
i
ci|φi〉|χi〉, (3.24)
which results in the reduced matrices
ρA =
i
|ci|^2 |φi〉〈φi|, ρB =
i
|ci|^2 |χi〉〈χi|. (3.25)
As the coefficients |ci|^2 are identical for both subsystems, the reduced density matrices have the same nonzero eigenvalues and therefore the same entropy. Therefore, we can define an “entanglement entropy” for any pure state |ψ〉 simply as the entropy of the reduced density matrix S(ρA). In general, the entanglement entropy will depend on the choice of the subsystems, however, we are mostly interested in the partitioning that maximizes the entanglement entropy. For translationally invariant systems, this is achieved simply by cutting the system in half. Coming back to the question on the usefulness of DMRG, we take a look at the entan- glement entropy of matrix product states. As the states are encoded in terms of a D × D matrix, the entanglement entropy of matrix product states satisfies the relation
S(ρA) ≤ 2 log D, (3.26)
