This page was generated from solutions_uniform_toycodes.ipynb (download).

Uniform matrix product states in the thermodynamic limit - Solutions

Closely following [1], we want to numerically investigate uniform matrix product states (uMPS) in the thermodynamic limit, importing the modules f_umps.py, g_vumps.py, h_utdvp.py and i_uexcitations.py.

Part 1) Injectivity and canonical form in f_umps.py

A uMPS \(\vert \psi(A) \rangle\) represents a translation invariant state of an infinite chain of \(d\)-dimensional quantum systems and is completely characterized by a single tensor \(A \in \mathbb{C}^{D \times d \times D}\) of physical dimension \(d\) and bond dimension \(D\):

\[\vert \psi(A) \rangle = \sum_{\{s\}} \langle v_L \vert \left(\prod_{n \in \mathbb{Z}} A^{s_n}\right) \vert v_R \rangle \vert \{s\} \rangle.\]

A crucial condition for efficient computations is injectivity of \(A\), meaning that counting algebraic multiplicities, the transfer matrix

\[T_A = \sum_{s=1}^d A^s \otimes \overline{A_s}\]

has only one eigenvalue of magnitude (and value) 1 and the corresponding right and left eigenvectors are positive definite matrices [2, 3]. Non-injective tensors appear with measure zero. In 1.1) we want to show that a randomly drawn tensor is indeed injective. Within the gauge freedom \(A^s \rightarrow XA^sX^{-1}\) for invertible \(X\), we can bring a uMPS in canonical form

\[\vert \psi(A) \rangle = \sum_{\{s\}} \cdot\cdot\cdot A_L^{s_{n-2}}A_L^{s_{n-1}}A_C^{s_n}A_R^{s_{n+1}}A_R^{s_{n+2}} \cdot\cdot\cdot \vert \{s\} \rangle,\]

with left/right orthonormal tensors \(A_L\)/\(A_R\) (injective with left/right leading eigenvector equal to identity matrix) and center site tensor \(A_C = A_L C = C A_R\) for some center matrix \(C\). These 4 tensors are the central attributes of the UniformMPS class. As a first physically meaningful and exactly representable uMPS, we consider the AKLT state and justify a few of its properties in 1.2). In 1.3) we make clear under which conditions a UniformMPS instance and a SimpleMPS(bc="infinite") instance (implemented in a_mps.py) can be converted into each other.

We then investigate the variational power of the set of uMPS with a given bond dimension \(D\), forming a manifold within the full Hilbert space. An element of the tangent space at point \(\vert \psi(A) \rangle\) in this manifold reads

\[\vert \psi(B; A) \rangle = \sum_n \sum_{\{s\}} \cdot\cdot\cdot A_L^{s_{n-2}}A_L^{s_{n-1}}B^{s_n}A_R^{s_{n+1}}A_R^{s_{n+2}} \cdot\cdot\cdot \vert \{s\} \rangle,\]

and we denote the projector onto the (orthogonal) tangent space with \(P_A\). For Hamiltonians \(H = \sum_n h_{n, n+1}\) with nearest-neighbor interaction \(h \in \mathbb{C}^{d^2 \times d^2}\), we run the following algorithms:

Part 2) Variational uniform matrix product states (VUMPS) in g_vumps.py

The tangent space gradient reads

\[P_A H \vert \psi(A) \rangle = \vert \psi(G; A)\rangle \:\:\mathrm{with}\:\: G = H_{\mathrm{eff},1}(A_C) - A_L H_{\mathrm{eff},0}(C) = H_{\mathrm{eff},1}(A_C) - H_{\mathrm{eff},0}(C) A_R.\]

The variational ground state optimum (corresponding to \(G=0\)) in the uMPS manifold satisfies a) \(A_C\) ground state of \(H_{\mathrm{eff},1}\), b) \(C\) ground state of \(H_{\mathrm{eff},0}\), c) \(A_C = A_L C = C A_R\). For this perform the following updates till convergence in gradient norm:

  1. \(A_C\) \(\rightarrow\) ground state of \(H_{\mathrm{eff},1}\), b) \(C\) \(\rightarrow\) ground state of \(H_{\mathrm{eff},0}\), c) \(A_L\)/\(A_R\) from left/right polar decompositions of \(A_C\) and \(C\). [1] >> Note that in contrast to iDMRG (implemented in d_dmrg.py), which successively grows the lattice by updated unit cells, VUMPS truly solves the variational problem in the sense of completely updating the state with each iteration and consequently keeping the translation invariance at any time. [4]

Part 3) Uniform time dependent variational principle (uTDVP) in h_utdvp.py

Approximately solves

\[\frac{d}{dt} \vert \psi(A) \rangle = -i P_A H \vert \psi(A) \rangle\]

for a small time step \(dt\) by a) \(A_C \rightarrow e^{-idtH_{\mathrm{eff},1}}(A_C)\), b) \(C \rightarrow e^{-idtH_{\mathrm{eff},0}}(C)\), c) \(A_L\)/\(A_R\) from left/right polar decompositions of \(A_C\) and \(C\). [1]

Part 4) Variational plane wave excitations in i_uexcitations.py

On top of a uMPS ground state \(\vert \psi(A) \rangle\), we want to variationally find quasiparticle excitations in a plane wave superposition of the form

\[\vert \psi(p,X; A) \rangle = \sum_n e^{ipn} \sum_{\{s\}} \cdot\cdot\cdot A_L^{s_{n-2}}A_L^{s_{n-1}}(V_L^{s_n}X)A_R^{s_{n+1}}A_R^{s_{n+2}} \cdot\cdot\cdot \vert \{s\} \rangle,\]

where \(p\) is the momentum, and \(X\) the left-gauge parametrization of the tensor \(B = V_L X\) (\(V_L\) is orthogonal to \(A_L\)), perturbing the ground state around site \(n\). With

\[\langle \psi(p',X'; A) \vert H \vert \psi(p,X; A) \rangle = 2\pi \delta(p'-p)\langle X' \vert H_{\mathrm{eff}}(p) \vert X \rangle,\]

the optimization boils down to diagonalizing \(H_{\mathrm{eff}}(p)\) for a few lowest-lying eigenvalues. [1, 5]

We benchmark our algorithms with the transverse field Ising (TFI) model, which can be diagonalized analytically in one dimension. [6]

References

[1] Vanderstraeten et al., Tangent-space methods for uniform matrix product states, 2019, https://arxiv.org/abs/1810.07006. [2] Perez-Garcia et al., Matrix product state representations, 2007, https://arxiv.org/abs/quant-ph/0608197. [3] Wolf, Quantum Channels and Operations - Guided Tour, 2012, https://mediatum.ub.tum.de/node?id=1701036. [4] Zauner-Stauber et al., Variational optimization algorithms for uniform matrix product states, 2018, https://arxiv.org/abs/1701.07035. [5] Haegeman et al., Variational matrix product ansatz for dispersion relations, 2012, https://arxiv.org/abs/1103.2286. [6] Subir Sachdev, Quantum Phase Transitions, 2nd ed, Cambridge University Press, 2011.

[1]:

import numpy as np
import matplotlib.pyplot as plt

1) f_umps.py: Random states and AKLT ground state

[2]:

from tenpy_toycodes.f_umps import UniformMPS, TransferMatrix
from tenpy_toycodes.a_mps import SimpleMPS

1.1) Transfer matrix and canonical form of random injective tensor

  • Create a random tensor \(A \in \mathbb{C}^{D \times d \times D}\) of physical dimension \(d = 2\) and bond dimension \(D = 4\).

[3]:

D, d = 4, 2
A = UniformMPS.get_random_tensor(D, d)

  • Show that \(A\) is injective, i.e. its TransferMatrix fulfills \(T_A = \vert R \rangle \langle L \vert + \mathcal{O}(\vert \lambda_2 \vert)\) with \(\vert \lambda_2 \vert < 1\) and \(R,L > 0\).

[4]:

# right eigenvector
T = TransferMatrix([A], [A])
lambdas, Rs = T.get_leading_eigenpairs(k=2)
lambda1, lambda2 = lambdas[0], lambdas[1]
R = Rs[0]
R /= np.trace(R)
ERs, _ = np.linalg.eig(R)

# left eigenvector
T = TransferMatrix([A], [A], transpose=True)
mus, Ls = T.get_leading_eigenpairs(k=2)
L = Ls[0]
L /= np.trace(L)
ELs, _ = np.linalg.eig(L)

assert np.allclose(np.abs(lambdas), np.abs(mus))
print("lambda_1 = ", lambda1)
print("|lambda_2| = ", np.abs(lambda2))
print("R > 0 with eigenvalues", np.real_if_close(ERs))
print("L > 0 with eigenvalues", np.real_if_close(ELs))

lambda_1 =  (0.9999999999999913-5.551115123125783e-17j)
|lambda_2| =  0.650921953008887
R > 0 with eigenvalues [0.50901127 0.30092656 0.05625866 0.13380351]
L > 0 with eigenvalues [0.63066192 0.22089672 0.03770395 0.11073741]

  • Compute the canonical form of \(A\), consisting of left/right orthonormal tensor \(A_L / A_R\), center site tensor \(A_C\), and center matrix \(C\).

[5]:

AL, AR, AC, C = UniformMPS.to_canonical_form(A)

AL, L: Converged up to tol=1e-10. Final error after 6 iterations: 2.6752155210835276e-13.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 9.133347930659786e-16.

  • Show that \(T_{A_L} = \vert R \rangle \langle \mathbb{1} \vert + \mathcal{O}(\vert \lambda_2 \vert)\) with \(\vert \lambda_2 \vert < 1\) and \(R > 0\).

[6]:

# right eigenvector
T = TransferMatrix([AL], [AL])
lambdas, Rs = T.get_leading_eigenpairs(k=2)
lambda1, lambda2 = lambdas[0], lambdas[1]
R = Rs[0]
R /= np.trace(R)
ERs, _ = np.linalg.eig(R)

# left eigenvector
T = TransferMatrix([AL], [AL], transpose=True)
mus, Ls = T.get_leading_eigenpairs(k=2)
L = Ls[0]
L /= np.trace(L)

assert np.allclose(np.abs(lambdas), np.abs(mus))
print("lambda_1 = ", lambda1)
print("|lambda_2| = ", np.abs(lambda2))
print("R > 0 with eigenvalues", np.real_if_close(ERs))
print("L = 1:", np.allclose(L, np.eye(4)/4))

lambda_1 =  (1.0000000000000038-1.6653345369377348e-16j)
|lambda_2| =  0.6509219530088591
R > 0 with eigenvalues [0.69095712 0.22999861 0.01121072 0.06783355]
L = 1: True

  • Show that \(T_{A_R} = \vert \mathbb{1} \rangle \langle L \vert + \mathcal{O}(\vert \lambda_2 \vert)\) with \(\vert \lambda_2 \vert < 1\) and \(L > 0\).

[7]:

# right eigenvector
T = TransferMatrix([AR], [AR])
lambdas, Rs = T.get_leading_eigenpairs(k=2)
lambda1, lambda2 = lambdas[0], lambdas[1]
R = Rs[0]
R /= np.trace(R)

# left eigenvector
T = TransferMatrix([AR], [AR], transpose=True)
lambdas, Ls = T.get_leading_eigenpairs(k=2)
lambda1, lambda2 = lambdas[0], lambdas[1]
L = Ls[0]
L /= np.trace(L)
ELs, _ = np.linalg.eig(L)

assert np.allclose(np.abs(lambdas), np.abs(mus))
print("lambda_1 = ", lambda1)
print("|lambda_2| = ", np.abs(lambda2))
print("R = 1:", np.allclose(R, np.eye(4)/4))
print("L > 0 with eigenvalues", np.real_if_close(ELs))

lambda_1 =  (1.000000000000003+0j)
|lambda_2| =  0.6509219530088584
R = 1: True
L > 0 with eigenvalues [0.69095712 0.22999861 0.01121072 0.06783355]

  • Show that \(A_C = A_L C = C A_R\).

[8]:

AC_L = np.tensordot(AL, C, axes=(2, 0))  # vL p [vR], [vL] vR -> vL p vR
AC_R = np.tensordot(C, AR, axes=(1, 0))  # vL [vR], [vL] p vR -> vL p vR
print(np.linalg.norm(AC_L - AC))
print(np.linalg.norm(AC_R - AC))

0.0
9.774997529390909e-16

  • Directly create a random UniformMPS (for different tolerances) and test its canonical form.

[9]:

psi_random = UniformMPS.from_desired_bond_dimension(D=100, d=3, tol=1.e-10)
print(psi_random.test_canonical_form())
psi_random = UniformMPS.from_desired_bond_dimension(D=100, d=3, tol=1.e-5)
print(psi_random.test_canonical_form())

AL, L: Converged up to tol=1e-10. Final error after 7 iterations: 5.160377003601066e-12.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 1.110476558538116e-14.
[2.84400437e-14 2.85188495e-14 0.00000000e+00 1.17873048e-14]
AL, L: Converged up to tol=1e-05. Final error after 6 iterations: 3.649522639547124e-06.
AR, R: Converged up to tol=1e-05. Final error after 6 iterations: 1.7342692594901902e-07.
[2.81704332e-14 2.90117149e-14 0.00000000e+00 1.73426926e-07]

1.2) AKLT state

For the spin-1 matrices \(S_x\), \(S_y\) and \(S_z\), we consider the two-site Hamiltonian matrix

\[h = S_x \otimes S_x + S_y \otimes S_y + S_z \otimes S_z + \frac{1}{3}(S_x \otimes S_x + S_y \otimes S_y + S_z \otimes S_z)^2 = 2P_{S=2}-2/3.\]

The Hamiltonian for multiple sites then simply reads \(H = \sum_{n} h_{n, n+1}\). The eigenenergy is minimized (with a value of \(-2/3\)), if every bond has total spin \(S = 0\) or \(S = 1\). This is fulfilled by distributing spin-1/2 singlets between all neighbor-sites and projecting every site onto the \(S = 1\) subspace. This is the AKLT state \(\vert \psi_{\mathrm{AKLT}} \rangle\). A few lines derivation shows that it can be represented exactly by a uMPS \(\vert \psi(A) \rangle\) with left orthonormal tensor

\[A_L^{+1} = \sqrt{\frac{2}{3}}\sigma^{+}, A_L^{0} = \sqrt{\frac{1}{3}}\sigma^{z}, A_L^{-1} = -\sqrt{\frac{2}{3}}\sigma^{-}.\]

  • Implement \(h, h^2 \in \mathbb{C}^{3 \times 3 \times 3 \times 3}\).

[10]:

d = 3

# spin-1 matrices
Sx = np.array([[0., 1., 0.],[1., 0., 1.],[0., 1., 0.]]) / np.sqrt(2)
Sy = -1.j * np.array([[0., 1., 0.],[-1., 0., 1.],[0., -1., 0.]]) / np.sqrt(2)
Sz = np.array([[1., 0., 0.],[0., 0., 0.],[0., 0., -1.]])

# Hamiltonian
S_dot_S = np.kron(Sx, Sx) + np.kron(Sy, Sy) + np.kron(Sz, Sz)
h = S_dot_S + (1/3.) * S_dot_S @ S_dot_S
h2 = h @ h
h = np.reshape(h, (d, d, d, d))
h2 = np.reshape(h2, (d, d, d, d))

  • Bring the above defined \(A_L\) into canonical form and initialize \(\vert \psi_{\mathrm{AKLT}} \rangle\) as a UniformMPS.

[11]:

D = 2

AL = np.zeros((D, d, D), dtype=complex)
AL[:, 0, :] = np.sqrt(2/3.) * np.array([[0., 0.],[1., 0.]])
AL[:, 1, :] = np.sqrt(1/3.) * np.array([[1., 0.],[0., -1.]])
AL[:, 2, :] = -np.sqrt(2/3.) * np.array([[0., 1.],[0., 0.]])

AR, C = UniformMPS.right_orthonormalize(AL)
AC = np.tensordot(AL, C, axes=(2, 0))  # vL p [vR], [vL] vR
psi_AKLT = UniformMPS(AL, AR, AC, C)

AR, R: Converged up to tol=1e-10. Final error after 6 iterations: 1.7071693894322978e-16.

Show the following properties of \(\vert \psi_{\mathrm{AKLT}} \rangle\):

  • Ground state energy \(e = \langle h \rangle = -2/3\),

  • Variance \(\langle h^2 \rangle - \langle h \rangle^2 = 0\) (H is frustration free),

  • Entanglement entropy \(S = \ln(2) \approx 0.6931471805599453\),

  • Correlation length \(\xi = -1/\ln(\vert \lambda_2 \vert) = -1/\ln(1/3) = 1/\ln 3 \approx 0.9102392266268371\),

  • Connected correlation function \(C(n) = \langle S_z^0 S_z^n \rangle - \langle S_z \rangle^2 \rightarrow \vert \lambda_2 \vert^{n-1}\).

[12]:

# ground state energy
print("e = ", psi_AKLT.get_bond_expectation_value(h))
# variance
print("var_e = ", psi_AKLT.get_bond_expectation_value(h2) - psi_AKLT.get_bond_expectation_value(h)**2)
# entanglement entropy
print("S = ", psi_AKLT.get_entanglement_entropy())
# correlation length
print("xi = ", psi_AKLT.get_correlation_length())
# second largest eigenvalue (to check consistency with formula above)
T = TransferMatrix([psi_AKLT.AL], [psi_AKLT.AL])
lambdas, _ = T.get_leading_eigenpairs(k=2)
lambda2 = lambdas[1]
print("|lambda2| = ", np.abs(lambda2))

e =  -0.6666666666666666
var_e =  -5.551115123125783e-17
S =  0.6931471805599453
xi =  0.9102392266268375
|lambda2| =  0.3333333333333333

[13]:

# correlation functions
N = 50
ns = np.arange(1, N+1)
Cs = psi_AKLT.get_correlation_functions(Sz, Sz, N)
Cs_connected = Cs - psi_AKLT.get_site_expectation_value(Sz)**2
Cs_theo = []
for n in ns:
    Cs_theo.append(np.abs(lambda2)**(n-1))

plt.figure(figsize=(6, 4))
plt.xlabel(r"$n$")
plt.semilogy(ns, Cs_connected, "x", label=r"$C(n)$")
plt.semilogy(ns, Cs_theo, "-", label=r"$\vert \lambda_2 \vert^{n-1}$")
plt.legend(loc="best")
plt.show()
../_images/toycodes_solutions_uniform_toycodes_27_0.png

1.3) Conversions between uniform MPS and infinite MPS

  • Convert the UniformMPS (uMPS) \(\vert \psi_{\mathrm{AKLT}} \rangle\) to a SimpleMPS(bc="infinite") (iMPS) and recheck the values of \(e\), \(S\), \(\xi\).

[14]:

psi_AKLT_infinite = psi_AKLT.to_infinite_MPS(L=1)

# ground state energy
print("e = ", psi_AKLT_infinite.bond_expectation_value([h]))
# variance
print("var_e = ", psi_AKLT_infinite.bond_expectation_value([h2]) \
      - psi_AKLT_infinite.bond_expectation_value([h])**2)
# entanglement entropy
print("S = ", psi_AKLT_infinite.entanglement_entropy())
# correlation length
print("xi = ", psi_AKLT_infinite.correlation_length())

e =  [-0.66666667]
var_e =  [-1.11022302e-16]
S =  [0.69314718]
xi =  0.9102392266268375

  • Denote by \(B \in \mathbb{C}^{2 \times 3 \times 2}\) the right canonical tensor of \(\vert \psi_{\mathrm{AKLT}} \rangle\) and by \(U \in \mathbb{C}^{2 \times 2}\) a random unitary. Show that the iMPS \(\vert\psi(B_1, B_2)\rangle\) with \(B_1 = BU\) and \(B_2 = U^*B\) is translation invariant and can be converted to a uMPS. Check that this is still the AKLT state.

[15]:

B = psi_AKLT_infinite.Bs[0]
C = psi_AKLT_infinite.Ss[0]
from scipy.stats import unitary_group
U = unitary_group.rvs(2)
B1 = np.tensordot(B, U, axes=(2, 0))
B2 = np.tensordot(np.conj(U).T, B, axes=(1, 0))

psi_infinite = SimpleMPS([B1, B2], [C, C], bc="infinite")
psi_uniform = UniformMPS.from_infinite_MPS(psi_infinite)

# ground state energy
print("e = ", psi_uniform.get_bond_expectation_value(h))
# variance
print("var_e = ", psi_uniform.get_bond_expectation_value(h2) - psi_uniform.get_bond_expectation_value(h)**2)
# entanglement entropy
print("S = ", psi_uniform.get_entanglement_entropy())
# correlation length
print("xi = ", psi_uniform.get_correlation_length())

Infinite MPS is translation invariant -> Conversion to uniform MPS.
AL, L: Converged up to tol=1e-10. Final error after 5 iterations: 6.61998662436351e-14.
AR, R: Converged up to tol=1e-10. Final error after 5 iterations: 6.072741058878367e-12.
e =  -0.666666666666667
var_e =  -3.3306690738754696e-16
S =  0.6931471805599452
xi =  0.9102392266268376

  • For random right canonical tensors \(B_1 \in \mathbb{C}^{D \times d \times D}\) and \(B_2 \in \mathbb{C}^{D \times d \times D}\), show that the iMPS \(\vert\psi(B_1, B_2)\rangle\) is in general not translation invariant and cannot be converted to a uMPS.

[16]:

tensor1 = UniformMPS.get_random_tensor(D=100, d=3)
tensor2 = UniformMPS.get_random_tensor(D=100, d=3)
_, B1, _, S1 = UniformMPS.to_canonical_form(tensor1)
_, B2, _, S2 = UniformMPS.to_canonical_form(tensor2)
S1, S2 = np.diag(S1), np.diag(S1)

psi_infinite = SimpleMPS([B1, B2], [S1, S2], bc="infinite")
psi_uniform = UniformMPS.from_infinite_MPS(psi_infinite)

AL, L: Converged up to tol=1e-10. Final error after 7 iterations: 1.2394265643359109e-11.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 1.1281732550398951e-14.
AL, L: Converged up to tol=1e-10. Final error after 7 iterations: 3.533688421854138e-12.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 1.3225822984078935e-14.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[16], line 8
      5 S1, S2 = np.diag(S1), np.diag(S1)
      7 psi_infinite = SimpleMPS([B1, B2], [S1, S2], bc="infinite")
----> 8 psi_uniform = UniformMPS.from_infinite_MPS(psi_infinite)

File ~/Library/Mobile Documents/com~apple~CloudDocs/master/thesis/python/tenpy_toycodes/tenpy_toycodes/f_umps.py:163, in UniformMPS.from_infinite_MPS(cls, psi)
    161 lambda1, U = T.get_leading_eigenpairs(k=1)
    162 if np.abs(np.abs(lambda1) - 1.) > 1.e-8:
--> 163     raise ValueError(f"Infinite MPS is not translation invariant.")
    164 else:
    165     print(f"Infinite MPS is translation invariant -> Conversion to uniform MPS.")

ValueError: Infinite MPS is not translation invariant.

Transverse field Ising model

In the following we want to use uMPS methods to find the ground state, elementary excitations and time-evolved states of the transverse field Ising (TFI) model

\[H = -J \sum_{n \in \mathbb{Z}} \sigma^z_n \sigma^z_{n+1} - g \sum_{n \in \mathbb{Z}} \sigma^x_n \overset{J=1}{=} -\sum_{n \in \mathbb{Z}} \sigma^z_n \sigma^z_{n+1} - g \sum_{n \in \mathbb{Z}} \sigma^x_n.\]

\(\mathbb{Z}_2-\mathrm{symmetry}\): \([H, U] = 0\) with \(U = \prod\limits_n \sigma^x_n\), local order parameter: magnetization density \(\langle \sigma^z \rangle\). From connecting the two limits of the transverse field \(g\) we expect the following quantum phase diagram:

  • Ferromagnetic phase \(g < g_c\)

    • two degenerate, symmetry broken ground states related by \(U\) (\(\ket{...\uparrow\uparrow\uparrow...}\) and \(\ket{...\downarrow\downarrow\downarrow...}\) for \(g = 0\))

    • \(\langle \sigma^z \rangle = \pm m \neq 0\) (\(m = 1\) for \(g = 0\))

    • elementary excitations: topological domain walls (\(\ket{...\uparrow\uparrow\uparrow \downarrow\downarrow\downarrow...}\) for \(g = 0\))

  • Paramagnetic phase \(g > g_c\)

    • unique symmetric ground state (\(\ket{...\rightarrow\rightarrow\rightarrow...}\) for \(g \to \infty\))

    • \(\langle \sigma^z \rangle = 0\)

    • elementary excitations: single spin flips (\(\ket{...\rightarrow\rightarrow \leftarrow \rightarrow\rightarrow...}\) for \(g \to \infty\))

  • Quantum phase transition at \(g_c\)

By performing Jordan-Wigner, Fourier and Bogoliubov transformations, the TFI model (with PBC) can be diagonalized analytically. The Hamiltonian in terms of fermionic creation and annihilation operators \(\gamma_{p}^{\dagger}\) and \(\gamma_{p}\) reads

\[H = \sum_{p} \epsilon_p \gamma_{p}^{\dagger}\gamma_{p} + E_0.\]

Single particle excitation energy: \(\epsilon_p = 2 \sqrt{1 - 2g\cos(p) + g^2}\) \(\Rightarrow\) energy gap closes at quantum critical point \(g_c = 1\).

Ground state energy: \(E_0 = -\sum_{p} \epsilon_p/2\) \(\Rightarrow\) ground state energy density in the thermodynamic limit: \(e_0 = - \frac{1}{2\pi} \int_{-\pi}^{\pi} \epsilon_p/2\). [2]

3) h_utdvp.py: Global quench dynamics

[22]:

from tenpy_toycodes.b_model import TFIModel
from tenpy_toycodes.c_tebd import example_TEBD_gs_tf_ising_infinite, calc_U_bonds, update_bond

from tenpy_toycodes.f_umps import UniformMPS
from tenpy_toycodes.g_vumps import vumps_algorithm
from tenpy_toycodes.h_utdvp import utdvp_algorithm

Consider the following global quench dynamics for the TFI model: the ground state for a value \(g_1\) of the transverse field is time-evolved according to a different value \(g_2\),

\[\ket{\psi(t)} = e^{-it H(g_2)} \ket{\psi_0(g_1)}.\]

We want to implement this time evolution using the utdvp_algorithm. For benchmark we use iTEBD. The run_TEBD function in c_tebd.py relies on a first order Trotter decomposition. We want to improve this to second order:

\[H = H_{\mathrm{odd}} + H_{\mathrm{even}} \Rightarrow U(dt) = e^{-i dt H_{\mathrm{odd}}/2} e^{-i dt H_{\mathrm{even}}} e^{-i dt H_{\mathrm{odd}}/2} + \mathcal{O}(dt^3).\]

  • Implement this improved scheme in a function run_TEBD_second_order.

  • Then write a function itebd_global_quench performing the global quench dynamics described above. Converge the method in bond dimension \(D\) and time step \(dt\).

[23]:

def run_TEBD_second_order(psi, U_bonds, U_bonds_half_dt, N_steps, chi_max, eps):
    """Evolve MPS psi for N_steps time steps with second order TEBD.

    H = H_odd + H_even
    U = exp(-i * dt * H) \\approx exp(-i * dt/2 * H_odd) exp(-i * dt * H_even) exp(-i* dt/2 *H_odd)
    """
    Nbonds = psi.L - 1 if psi.bc == 'finite' else psi.L
    assert len(U_bonds) == len(U_bonds_half_dt) == Nbonds
    for i_bond in range(0, Nbonds, 2):  # odd (half)
        update_bond(psi, i_bond, U_bonds_half_dt[i_bond], chi_max, eps)
    for i_bond in range(1, Nbonds, 2):  # even
        update_bond(psi, i_bond, U_bonds[i_bond], chi_max, eps)
    for n in range(N_steps - 1):  # N_steps - 1 times odd, even
        for k in [0, 1]:
            for i_bond in range(k, Nbonds, 2):
                update_bond(psi, i_bond, U_bonds[i_bond], chi_max,eps)
    for i_bond in range(0, Nbonds, 2):  # odd (half)
        update_bond(psi, i_bond, U_bonds_half_dt[i_bond], chi_max, eps)

def itebd_global_quench(g1, g2, chi_max, dt, T):
    """Calculate the tfi ground state for transverse field g1 and evolve it according to different g2
    up to time T in steps of dt with second order iTEBD."""
    _, psi1, _ = example_TEBD_gs_tf_ising_infinite(g1, chi_max)
    model2 = TFIModel(L=2, J=1., g=g2, bc="infinite")
    H_bonds = model2.H_bonds
    U_bonds = calc_U_bonds(H_bonds, 1.j * dt)
    U_bonds_half_dt = calc_U_bonds(H_bonds, 1.j * dt/2)
    ts = []
    Ss = []
    t = 0.
    ts.append(t)
    Ss.append(psi1.entanglement_entropy()[1])
    psi = psi1
    for i in range(int(T/dt)):
        run_TEBD_second_order(psi, U_bonds, U_bonds_half_dt, N_steps=1, chi_max=chi_max, eps=1.e-10)
        t += dt
        ts.append(t)
        Ss.append(psi.entanglement_entropy()[0])
    return ts, Ss

[24]:

g1 = 3.
g2 = 1.5
T = 5.

[25]:

plt.figure(figsize=(6, 4))
plt.title(f"TFI global quench with $g_2={g2}$ on GS of $g_1={g1}$ (iTEBD)")
plt.xlabel('time $t$')
plt.ylabel('entanglement entropy')
for dt in [0.1, 0.01, 0.001]:
    for D in [10, 20, 30]:
        ts_itebd, Ss_itebd = itebd_global_quench(g1, g2, D, dt, T)
        plt.plot(ts_itebd, Ss_itebd, label=f"D={D}, dt={dt}")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()

infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [10, 10]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.886591626230449
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [11, 11]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8939157106447442
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [11, 11]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8939157106447435
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [10, 10]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8865916262304465
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [11, 11]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8939157106447669
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [11, 11]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8939157106447387
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [10, 10]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8865916262304483
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [11, 11]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8939157106447684
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [11, 11]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8939157106447557
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
../_images/toycodes_solutions_uniform_toycodes_49_1.png

  • Analogous to itebd_global_quench, write a function utdvp_global_quench and make sure they give the same dynamics.

[26]:

def utdvp_global_quench(g1, g2, D, dt, T):
    """Calculate the tfi ground state for transverse field g1 and evolve it according to different g2
    up to time T in steps of dt with uTDVP."""
    h1 = TFIModel(L=1, J=1., g=g1, bc="infinite").H_bonds[0]
    h2 = TFIModel(L=1, J=1., g=g2, bc="infinite").H_bonds[0]
    guess_psi1 = UniformMPS.from_desired_bond_dimension(D)
    _, psi1, _ = vumps_algorithm(h1, guess_psi1, tol=1.e-10)
    ts, Ss = utdvp_algorithm(psi1, h2, dt, T)
    return ts, Ss

[27]:

g1 = 3.
g2 = 1.5
D = 30
T = 5.
dt = 0.01

[28]:

ts_itebd, Ss_itebd = itebd_global_quench(g1, g2, D, dt, T)
ts_utdvp, Ss_utdvp = utdvp_global_quench(g1, g2, D, dt, T)

infinite TEBD, imaginary time evolution, transverse field Ising
g=3.00
dt = 0.10000: E (per site) = -3.0722915530137
dt = 0.01000: E (per site) = -3.0827483842174
dt = 0.00100: E (per site) = -3.0838108901537
dt = 0.00010: E (per site) = -3.0839168389419
dt = 0.00001: E (per site) = -3.0839274447472
final bond dimensions:  [11, 11]
<sigma_x> = 0.97161
<sigma_z> = -0.00000
correlation length: 0.8939157106447625
Analytic result: E (per site) = -3.0839288503801
relative error:  4.5579290330811456e-07
AL, L: Converged up to tol=1e-10. Final error after 7 iterations: 4.989161815670688e-11.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 3.2188988713804404e-15.
uMPS ground state converged with VUMPS up to tol=1e-10 in gradient norm. Final error after 4 iterations: 2.4714217379505394e-12.
Ground state energy density: -3.083928850380077.
Ground state variance density: (-6.245004513516506e-16-5.5294310796760726e-17j).
uMPS evolved with TDVP up to time T=5.0 in steps of dt=0.01.

[29]:

plt.figure(figsize=(6, 4))
plt.title(f"TFI global quench with $g_2={g2}$ on GS of $g_1={g1}$")
plt.xlabel('time $t$')
plt.ylabel('entanglement entropy')
plt.plot(ts_itebd, Ss_itebd, label=f"iTEBD (D={D}, dt={dt})")
plt.plot(ts_utdvp, Ss_utdvp, label=f"uTDVP (D={D}, dt={dt})")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()
../_images/toycodes_solutions_uniform_toycodes_54_0.png

4) i_uexcitations.py: Variational plane wave excitations

[30]:

from tenpy_toycodes.tfi_exact import infinite_gs_energy, infinite_excitation_dispersion
from tenpy_toycodes.b_model import TFIModel

from tenpy_toycodes.f_umps import UniformMPS
from tenpy_toycodes.g_vumps import vumps_algorithm
from tenpy_toycodes.i_uexcitations import VariationalPlaneWaveExcitationEngine

By running the VariationalPlaneWaveExcitationEngine, we want to compute the single particle excitations \(\epsilon_p = 2 \sqrt{1 - 2g\cos(p) + g^2}\) on top of the TFI ground state.

  • For this write a function get_tfi_spectrum, returning both the variational dispersion relation and the exact one from tfi_exact.py. Actively target topological domain wall excitations for transverse field values \(g < 1\).

[31]:

def get_tfi_spectrum(g, D):
    print(f"g = {g:.1f}, D = {D}:")
    # exact analytical solution
    e0_exact = infinite_gs_energy(J=1., g=g)
    print(f"e0_exact = {e0_exact}.")
    ps_exact, es_exact = infinite_excitation_dispersion(J=1., g=g)
    # vumps and variational plane wave excitations
    h = TFIModel(L=1, J=1., g=g, bc='infinite').H_bonds[0]
    guess_psi0 = UniformMPS.from_desired_bond_dimension(D)
    e0, psi0, _ = vumps_algorithm(h, guess_psi0, tol=1.e-8)
    print(f"|e0_exact - e0| = {np.abs(e0_exact - e0)}")
    ps = list(np.arange(-np.pi, 0, np.pi/10)) + [0] + list(np.arange(0, np.pi + np.pi/10, np.pi/10))
    if g >= 1.:  # paramagnetic phase (single spin flip excitations) and critical point (gapless excitation)
        es = []
        for p in ps:
            excitation_engine = VariationalPlaneWaveExcitationEngine(psi0, h, p, k=1, tol=1.e-8)
            e, _ = excitation_engine.run()
            es.append(e)
        return ps_exact, es_exact, ps, es
    elif 0 <= g < 1:  # ferromagnetic phase (topological domain wall excitations)
        sigma_z = TFIModel(L=1, J=1., g=g, bc='infinite').sigmaz
        sz = psi0.get_site_expectation_value(sigma_z)
        for i in range(100):
            guess_psi0_tilde = UniformMPS.from_desired_bond_dimension(D)
            _, psi0_tilde, _ = vumps_algorithm(h, guess_psi0_tilde, tol=1.e-8)
            sz_tilde = psi0_tilde.get_site_expectation_value(sigma_z)
            if np.abs(sz + sz_tilde) < 1.e-5:
                print(f"Found sz={sz} and sz_tilde={sz_tilde} after {i+2} ground state searchs.")
                es = []
                for p in ps:
                    excitation_engine = VariationalPlaneWaveExcitationEngine(psi0, h, p, k=1, tol=1.e-8, \
                                                                             psi0_tilde=psi0_tilde)
                    e, _ = excitation_engine.run()
                    es.append(e)
                return ps_exact, es_exact, ps, es

[32]:

D = 30
ps_exact, es_exact_para, ps, es_para = get_tfi_spectrum(g=1.5, D=D)
_, es_exact_crit, _, es_crit = get_tfi_spectrum(g=1.0, D=D)
_, es_exact_ferro, _, es_ferro = get_tfi_spectrum(g=0.5, D=D)

g = 1.5, D = 30:
e0_exact = -1.6719262215361952.
AL, L: Converged up to tol=1e-10. Final error after 7 iterations: 3.960938757655051e-11.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 6.288235903606723e-14.
uMPS ground state converged with VUMPS up to tol=1e-08 in gradient norm. Final error after 4 iterations: 3.897649278681079e-09.
Ground state energy density: -1.671926221535125.
Ground state variance density: (-3.0839913955915677e-12+4.542807102714264e-16j).
|e0_exact - e0| = 1.070254995738651e-12
g = 1.0, D = 30:
e0_exact = -1.2732395447351625.
AL, L: Converged up to tol=1e-10. Final error after 7 iterations: 8.370035562611671e-12.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 1.3174723448174354e-15.
uMPS ground state converged with VUMPS up to tol=1e-08 in gradient norm. Final error after 78 iterations: 9.892290565137428e-09.
Ground state energy density: -1.2732395297253791.
Ground state variance density: (2.3240924379441807e-08+1.6151901864602802e-16j).
|e0_exact - e0| = 1.5009783416175537e-08
g = 0.5, D = 30:
e0_exact = -1.0635444099733717.
AL, L: Converged up to tol=1e-10. Final error after 7 iterations: 7.243643498737732e-12.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 1.3140283050486809e-15.
uMPS ground state converged with VUMPS up to tol=1e-08 in gradient norm. Final error after 4 iterations: 1.6905444333507086e-11.
Ground state energy density: -1.063544409973369.
Ground state variance density: (1.8409752888803865e-16-1.0965688535154172e-17j).
|e0_exact - e0| = 2.6645352591003757e-15
AL, L: Converged up to tol=1e-10. Final error after 7 iterations: 3.590882800043658e-13.
AR, R: Converged up to tol=1e-10. Final error after 7 iterations: 1.266832377228487e-15.
uMPS ground state converged with VUMPS up to tol=1e-08 in gradient norm. Final error after 3 iterations: 5.003229193865994e-09.
Ground state energy density: -1.0635444099732108.
Ground state variance density: (-2.7150373962947505e-13+7.430173013314723e-18j).
Found sz=0.9646786299603096 and sz_tilde=-0.9646786299571539 after 2 ground state searchs.

[33]:

plt.figure(figsize=(6, 4))
plt.title(r"TFI excitations $\epsilon_p = 2 \sqrt{g^2 - 2g\cos(p) + 1}$")
plt.xticks(np.arange(-np.pi, 3*np.pi/2, np.pi/2), \
           [r"$-\pi$",r"$-\pi /2$", r"$0$", r"$\pi /2$",r"$\pi$"])
plt.axhline(y=0., color="black")
plt.plot(ps_exact, es_exact_para, color="lightskyblue")
plt.plot(ps, es_para, ".", color="blue", label=r"$g = 1.5$")
plt.plot(ps_exact, es_exact_crit, color="lightcoral")
plt.plot(ps, es_crit, ".", color="darkred", label=r"$g = 1.0$")
plt.plot(ps_exact, es_exact_ferro, color="lightgreen")
plt.plot(ps, es_ferro, ".", color="green", label=r"$g = 0.5$")
plt.xlabel(r"$p$")
plt.ylabel(r"$\epsilon_p$")
plt.legend(title=f"VUMPS + variational perturbations ($D = {D}$)", loc='center left', \
           bbox_to_anchor=(1, 0.5))
plt.show()
../_images/toycodes_solutions_uniform_toycodes_60_0.png