Probing Many-body Localization by Excited-state VQE#

Overview#

This tutorial introduces a new method to probe many-body localization (MBL) by excited-state VQE. The model hosting MBL transition considered here is the interacting Aubry-Andr\(\acute{e}\) model which reads:

\[\begin{split}H=\sum_{i}(\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1} +V_{0}\sigma_{i}^{z} \sigma_{i+1}^{z} ) \\ +W_{0} \sum_{i=1}^L \cos(2 \pi \eta i+\phi) \sigma_{i}^{z},\end{split}\]

where \(\sigma^{x,y,z}_{i}\) are Pauli matrices of the \(i\)-th qubit, \(L\) is the number of total qubits, \(V_{0}\) is the \(zz\) interaction strength, \(W\) is the strength of quasi-periodic (QP) potential and \(\phi\) is the phase of the cosine potential. We set \(V_{0}=0.5\) and \(\eta=(\sqrt{5}-1)/2\). When \(W\) is large, the system enters the MBL phase and exhibits exotic behaviors, such as “area law” entanglement for highly excited states.

Setup#

[1]:

import time
import numpy as np
import tensorflow as tf
import tensorcircuit as tc

You can define different gates easily in TensorCircuit.

[2]:

tc.set_backend("tensorflow")
tc.set_dtype("complex128")
dtype = np.complex128

ii = tc.gates._ii_matrix
xx = tc.gates._xx_matrix
yy = tc.gates._yy_matrix
zz = tc.gates._zz_matrix

Note: The implementation in this tutorial is fixed to the TensorFlow backend. One can use the APIs of TensorCircuit and the APIs of TensorFlow seamlessly and freely, and the automatic differential and the JIT compilation will not be influenced.

Parameters#

[3]:

L = 4
V0 = 0.5
W0 = [1.0, 10.0]  # W0=1.0: thermal phase; W0=10.0: MBL phase
η = (np.sqrt(5) - 1) / 2
ϕ = 0.1

num_trial = 5  # The number of independent trials
depth = 1  # The depth of the PQC
θ0 = np.pi / 5.0  # The rotation angle of the input state

Construct Circuit#

The circuit structure for the excited-state VQE and eigenstate witness measurement. Here \(U_{0}(\theta)= \rm{exp}(\frac{i\theta}{2}(\sigma^{x} \otimes \sigma^{x} + \sigma^{y} \otimes \sigma^{y}))\).

total.jpeg

[4]:

@tf.function
def get_ρ_circuit(v, epochs=2):
    c = tc.DMCircuit(L)

    ### Part I: the input-state preparation quantum circuit
    ###         prepare an antiferromagnetic (AF) state (|0101...>) from the initial state |0000...>
    ###         and the system will stay in Mz=0 sector due to U(1) symmetry
    for i in range(L):
        if i % 2 == 1:
            c.ry(i, theta=tf.cast(np.pi, dtype=dtype))

    ###         prepare input state
    for i in range(L):
        c.any(
            i,
            (i + 1) % L,
            unitary=0.5 * tf.math.cos(tf.cast(θ0, dtype=dtype)) * (ii - zz)
            + 0.5 * (ii + zz)
            + 0.5 * tc.backend.i() * tf.math.sin(tf.cast(θ0, dtype=dtype)) * (xx + yy),
        )  ### U_{0}(θ)
    ### Part II: the parameterized quantum circuit (PQC): optimize the parameters based on gradient descent
    for epoch in range(epochs):
        for i in range(L):
            c.any(
                i,
                (i + 1) % L,
                unitary=0.5
                * tf.math.cos(tf.cast(v[2 * epoch, i], dtype=dtype))
                * (ii - zz)
                + 0.5 * (ii + zz)
                + 0.5
                * tc.backend.i()
                * tf.math.sin(tf.cast(v[2 * epoch, i], dtype=dtype))
                * (xx + yy),
            )  ### U_{0}(θ)

        for i in range(L):
            c.rz(i, theta=tf.cast(v[2 * epoch + 1, i], dtype=dtype))
    ρ = c.densitymatrix()
    return ρ

Construct Hamiltonian#

[5]:

def get_hamiltonian(W0):
    h = []
    w = []
    ###XX
    for i in range(L):
        h.append([])
        w.append(1.0)
        for j in range(L):
            if j == (i + 1) % L or j == i:
                h[i].append(1)
            else:
                h[i].append(0)
    ###YY
    for i in range(L):
        h.append([])
        w.append(1.0)
        for j in range(L):
            if j == (i + 1) % L or j == i:
                h[i + L].append(2)
            else:
                h[i + L].append(0)
    ###ZZ
    for i in range(L):
        h.append([])
        w.append(V0)
        for j in range(L):
            if j == (i + 1) % L or j == i:
                h[i + 2 * L].append(3)
            else:
                h[i + 2 * L].append(0)
    ###potential
    for i in range(L):
        h.append([])
        w.append(W0 * np.cos(2 * np.pi * η * i + ϕ))
        for j in range(L):
            if j == i:
                h[i + 3 * L].append(3)
            else:
                h[i + 3 * L].append(0)

    hamiltonian = tc.quantum.PauliStringSum2Dense(
        tf.constant(h, dtype=tf.complex128), tf.constant(w, dtype=tf.complex128)
    )
    return hamiltonian

Loss Function#

The loss function used here is the energy variance that vanishes for eigenstates:

\[C(\theta) = \langle H^{2} \rangle - \langle H \rangle ^2\]
[6]:

@tf.function
def variance_loss(hamiltonian, v, epochs=2):
    with tf.GradientTape() as tape:
        tape.watch(v)
        ρ = get_ρ_circuit(v, epochs)
        hexp = tf.linalg.trace(ρ @ hamiltonian)
        loss = tf.linalg.trace(ρ @ hamiltonian @ hamiltonian) - hexp**2
        loss = tf.math.real(loss)
    gr = tape.gradient(loss, v)
    return loss, gr, ρ

Eigenstate Witness#

Under the excited-state VQE, the output states will converge to some highly excited states, and the MBL phase has a higher eigenstate witness.

[7]:

@tf.function
def get_eigenstate_witness(ρ, hamiltonian, evolution_time):
    evolution_p = tf.linalg.expm(
        tc.backend.i() * tf.cast(hamiltonian, dtype=dtype) * evolution_time
    )
    evolution_m = tf.linalg.expm(
        tc.backend.i() * tf.cast(hamiltonian, dtype=dtype) * (-evolution_time)
    )

    ρ_00 = tf.linalg.trace(ρ)
    ρ_01 = tf.linalg.trace(ρ @ evolution_p)
    ρ_10 = tf.linalg.trace(evolution_m @ ρ)
    ρ_11 = tf.linalg.trace(evolution_m @ ρ @ evolution_p)

    ρ_cq_up = tf.stack([ρ_00, ρ_01], axis=0)
    ρ_cq_dn = tf.stack([ρ_10, ρ_11], axis=0)
    ρ_cq = 0.5 * tf.stack([ρ_cq_up, ρ_cq_dn], axis=1)
    r = tf.linalg.trace(ρ_cq @ ρ_cq)
    return tf.math.real(r)

[8]:

def opt_main(hamiltonian, epochs=2, maxiter=200):
    v = tf.Variable(
        initial_value=tf.random.normal(
            shape=[2 * epochs, L], mean=0.0, stddev=0.03, dtype=tf.float64
        )
    )  # initialize parameters
    opt = tf.keras.optimizers.Adam(learning_rate=0.02)
    print("loop begins")
    for i in range(maxiter):
        loss, gr, ρ = variance_loss(hamiltonian, v.value(), epochs)
        opt.apply_gradients([(tf.math.real(gr), v)])
        if i == maxiter - 1:
            print(loss.numpy())
    return ρ

[9]:

start = time.time()
ρ_list = []
r_list = []
for i in range(len(W0)):
    hamiltonian = get_hamiltonian(W0[i])
    evolution_time = 10.0 / W0[i]
    ρ_list.append([])
    r_list.append([])
    for j in range(num_trial):
        ρ = opt_main(hamiltonian, depth)
        ρ_list[i].append(ρ)
        r = get_eigenstate_witness(ρ, hamiltonian, evolution_time)
        r_list[i].append(r)
endtime = time.time()
print("Run time = ", endtime - start)
print(
    "Eigenstate Witness:",
    np.mean(r_list[0]),
    "(Thermal Phase);",
    np.mean(r_list[1]),
    "(MBL Phase)",
)

loop begins
0.18502121452649556
loop begins
0.3647578794270281
loop begins
1.1563841018010805
loop begins
1.1626719640361998
loop begins
0.20636721189392704
loop begins
0.40186462162307635
loop begins
0.39881574091884886
loop begins
0.41128232341327475
loop begins
0.4055560658837294
loop begins
0.4046518612466343
Run time =  16.519852876663208
Eigenstate Witness: 0.8140801354126028 (Thermal Phase); 0.9994227942365521 (MBL Phase)

Research Project#

Reference paper: https://arxiv.org/pdf/2111.13719.pdf.