Quimb circuit based Trotterization vs TEBD in-built - result mismatch

[sajjan02purdue]

Hi my issue is the following.

System setup
################
I am trying to measure the following quantity
$F(t) = |\langle \psi |e^{-iHt} |\psi\rangle|^2$

where the hamiltonian $H$ is 1D XXZ chain and is defined as a HamLocal1D object generated by calling Hamiltonian_generation() function defined below.

1. $|\psi\rangle$ is a randomly prepared state from a quimb circuit generated by calling the function random_state_prep_circuit_quimb() defined below.

I am implementing $e^{-i H t}$ by two different means .

2. Using in-built TEBD (by calling the function in_built_TEBD() defined below) which is internal to quimb. I take the MPS from TEBD using at_times() command and then contract it against the initial random MPS (see above in point 1).

3. Using a circuit-based construction of second-order Trotterization of 1D XXZ obtained from a paper (snapshot shared below) . The circuit is correct and has been replicated in many papers.

I generate the circuit object in quimb , apply relevant gates and then extract the circuit MPS and then contract it against the initial random MPS (see above in point 1). I am using just 4 qubits for testing now and the parameters of XXZ is Jx=Jy=Jz=J. I have set dt = 0.1/J upto 40/J total timesteps . J is set to 1.

Issue/Problem statement
####################

While 2) above works fine and gives expected outcomes, the fidelities $F(t)$ generated by 3) are way too low and decreases extremely fast in time. I would like to know what is the reason for the discrepancy between the two methods
. I have tried many different things including increasing the bond dimension of the MPS representation of the circuit but to no avail. I would just like to ensure if I am doing something odd (which I suspect). All code required for minimal reproduction of the issue is provided below

Any help will be greatly appreciated.

Code and Function base

#####################

def Hamiltonian_generation(system_size, spin_type=1/2, bound_cond='obc', coeff_vector=[1,1,1], 
                           on_site=None):
    L=system_size
    builder = qtn.SpinHam1D(S=spin_type)
        
    if bound_cond == 'obc':
        connectivity_list = [(i,i+1) for i in np.arange(L-1)]
    else:
        connectivity_list = [(i,i+1) for i in np.arange(L-1)] + [(L-1, 0)]
        
    for qubit_tuple in connectivity_list:
        builder[qubit_tuple[0], qubit_tuple[1]] += coeff_vector[0], 'X', 'X'
        builder[qubit_tuple[0], qubit_tuple[1]] += coeff_vector[1], 'Y', 'Y'
        builder[qubit_tuple[0], qubit_tuple[1]] += coeff_vector[2], 'Z', 'Z'
        
    if on_site != None:
        for qubit in np.arange(L):
            builder[qubit] += on_site[0], 'X'
            builder[qubit] += on_site[1], 'Y'
            builder[qubit] += on_site[2], 'Z'

    H = builder.build_local_ham(L)
    return H


from scipy.stats import rv_continuous

class sin_sampler_dist(rv_continuous):
    def _pdf(self, theta):
        return 0.5*np.sin(theta)
sin_sampler = sin_sampler_dist(a=0, b=np.pi)


def random_state_prep_circuit_quimb(circuit, system_size, no_layers=3):
    from numpy.random import uniform
    #circ = qtn.Circuit(system_size, tag_gate_numbers=True, tag_gate_rounds=True, tag_gate_labels=True)
    regs = list(range(system_size))
    sin_sampler = sin_sampler_dist(a=0, b=np.pi)
    even_odd_pair_list = [[(i,i+1) for i in np.arange(0, system_size-1, 2)], [(i,i+1) for i in np.arange(1, system_size-1, 2)]]
    for _ in range(no_layers):
        for connectivity_list in even_odd_pair_list:
            for pair in connectivity_list:
                circuit.apply_gate('RZ', uniform(-np.pi, np.pi), regs[pair[0]])
                circuit.apply_gate('RY', float(sin_sampler.rvs(size=1)), regs[pair[0]])
                circuit.apply_gate('RZ', uniform(-np.pi, np.pi), regs[pair[0]])            

                circuit.apply_gate('RZ', uniform(-np.pi, np.pi), regs[pair[1]])
                circuit.apply_gate('RY', float(sin_sampler.rvs(size=1)), regs[pair[1]])
                circuit.apply_gate('RZ', uniform(-np.pi, np.pi), regs[pair[1]])
                
                circuit.apply_gate('CZ', regs[pair[0]], regs[pair[1]])  
            #circ.apply_gate('CNOT', regs[pair[0]], regs[pair[1]])
    return circuit


def unitary_2_qubit_infinitesimal_quimb(circuit, dt, coeff_vector, pair):
    import quimb
    theta_x = dt*coeff_vector[0]*(0.25**0)
    theta_y = dt*coeff_vector[1]*(0.25**0)
    theta_z = dt*coeff_vector[2]*(0.25**0)

    circuit.apply_gate('CX',pair[0],pair[1])
    circuit.apply_gate('RX',2*theta_x, pair[0])
    circuit.apply_gate('RZ',2*theta_z, pair[1])

    circuit.apply_gate('H',pair[0])
    circuit.apply_gate('CX',pair[0], pair[1])
    #circuit.apply_gate('RZ',-np.pi/2, pair[0])
    circuit.apply_gate("PHASE", -np.pi/2, pair[0])
    
    circuit.apply_gate('H',pair[0])
    circuit.apply_gate('RZ',-2*theta_y, pair[1])

    circuit.apply_gate('CX',pair[0], pair[1])
    circuit.apply_gate('RX',np.pi/2, pair[0])
    circuit.apply_gate('RX',-np.pi/2, pair[1])

    return circuit
    #return quimb.core.dag(circuit.uni)



def time_evolution_finite_quimb(quimb_circ, system_size, dt, No_trotter_steps, coeff_vector, on_site=[0,0,0]):
    
    for step_index in np.arange(1, No_trotter_steps+1):
        if on_site not in [None, [0,0,0]]:
            for qubit in np.arange(0, system_size):
                unitary_1_qubit_infinitesimal_quimb(quimb_circ, qubit, delta_t=dt/2.0, on_site=on_site)
                
        for pair in [(i,i+1) for i in np.arange(0, system_size, 2)]:
            quimb_circ = unitary_2_qubit_infinitesimal_quimb(quimb_circ, dt/2.0, coeff_vector, pair)

        for pair in [(i,i+1) for i in np.arange(1, system_size-1, 2)]:
            quimb_circ = unitary_2_qubit_infinitesimal_quimb(quimb_circ, dt, coeff_vector, pair)

        for pair in [(i,i+1) for i in np.arange(0, system_size, 2)]:
            quimb_circ = unitary_2_qubit_infinitesimal_quimb(quimb_circ, dt/2.0, coeff_vector, pair)

        if on_site not in [None, [0,0,0]]:
            for qubit in np.arange(0, system_size):
                unitary_1_qubit_infinitesimal_quimb(quimb_circ, qubit, delta_t=dt/2.0, on_site=on_site)
            
    return quimb_circ


def in_built_TEBD(Ham_Local_Ham1D, initial_state_MPS, 
                  total_time_list, dt=1e-2):
    
    tebd = qtn.TEBD(initial_state_MPS, Ham_Local_Ham1D, dt=dt, split_opts=dict(cutoff=1e-8), imag=False)

    #ts = np.linspace(0.0, 10.0, 51)
    fidelity_dict = {}
    fidelity_dict[0] = (abs(initial_state_MPS.H @ initial_state_MPS)**2)

    for index, psi_t in enumerate(tebd.at_times(total_time_list, dt=dt, order=2, progbar=False)):
        
        fidelity_dict[total_time_list[index]] = (abs(initial_state_MPS.H @ psi_t)**2)

    return fidelity_dict

Results

########

#System attributes and Hamiltonian generation (as LocalHam1D object)
###################################
system_size=4
spin_type=1/2
on_site=[0.0,0.0,0.0]
coeff_vector = [1,1,1]
no_random_states = 1
max_bond_dim = 20
Hamiltonian_params = {'bound_cond' : 'obc', 'coeff_vec' : coeff_vector, 'on-site' : on_site}
total_no_steps = 40
dt = 0.1
inds = 'd_{}'


Ham_Local_Ham1D = Hamiltonian_generation(system_size=system_size, spin_type=1/2, 
                           bound_cond=Hamiltonian_params['bound_cond'], 
                           coeff_vector=Hamiltonian_params['coeff_vec'], 
                           on_site=Hamiltonian_params['on-site'])



#Initial random state prep
#*************************
circuit = qtn.Circuit(system_size)
regs = list(range(system_size))
        
quimb_rand_circuit = random_state_prep_circuit_quimb(circuit, system_size=system_size, no_layers=1)
quimb_rand_circuit_MPS = qtn.CircuitMPS.from_gates(gates=quimb_rand_circuit.gates,
                              gate_opts=dict(
                              max_bond=10*max_bond_dim,
                              cutoff=1e-8),
                              progbar=False,
                               ).psi

#In-built TEBD call with initial random state
#***********************************
ovlp_inbuilt = in_built_TEBD(Ham_Local_Ham1D=Ham_Local_Ham1D, initial_state_MPS=quimb_rand_circuit_MPS, 
                        total_time_list=np.round(np.multiply(dt, np.arange(1, total_no_steps+1)), 3), dt=dt)
print ("ovlp inbuilt", ovlp_inbuilt)


#quimb circuit TEBD call with initial random state
#************************************************
for time in tqdm(np.round(np.multiply(dt, np.arange(1, total_no_steps+1)), 3)):
    quimb_rand_circuit = time_evolution_finite_quimb(quimb_rand_circuit, system_size, dt, No_trotter_steps=1, coeff_vector = coeff_vector, on_site=on_site)
    ovlp_quimb_circ[time] = np.abs(quimb_rand_circuit.psi.H @ quimb_rand_circuit_MPS)**2

print ("ovlp quimb circuit", ovlp_quimb_circ)


Obtained values are the following :

ovlp inbuilt 
{0: 0.9999999999999996, 0.1: 0.9959420262043145, 0.2: 0.983891407115981, 0.3: 0.9642130056409608, 0.4: 0.9374983467005762, 0.5: 0.9045415526145574, 0.6: 0.8663072221496924, 0.7: 0.8238917291899889, 0.8: 0.778479710461427, 0.9: 0.7312977182400782, 1.0: 0.6835671237596156, 1.1: 0.6364583653384801, 1.2: 0.5910485424703097, 1.3: 0.5482841687956572, 1.4: 0.5089506233444718, 1.5: 0.4736494953004309, 1.6: 0.44278462084450815, 1.7: 0.4165571819382002, 1.8: 0.39496979816803934, 1.9: 0.37783911620853705, 2.0: 0.36481600838146555, 2.1: 0.3554121514988259, 2.2: 0.34903148598357875, 2.3: 0.3450048656601664, 2.4: 0.3426261086650983, 2.5: 0.34118765292126413, 2.6: 0.34001410391814757, 2.7: 0.33849213176064225, 2.8: 0.33609541785946806, 2.9: 0.33240365478653927, 3.0: 0.3271149483662977, 3.1: 0.3200513397526423, 3.2: 0.31115753692245507, 3.3: 0.3004932998064056, 3.4: 0.2882202425838816, 3.5: 0.2745840841878343, 3.6: 0.25989358062978857, 3.7: 0.24449750100147927, 3.8: 0.22876105782765893, 3.9: 0.2130431711909683, 4.0: 0.19767583851098744}

100%|███████████████████████████████████████████| 40[/40]

ovlp quimb circuit 
{0: 0.9999999999997653, 0.1: 0.9369746964571253, 0.2: 0.7765595355083472, 0.3: 0.5873344146119812, 0.4: 0.43750202815996575, 0.5: 0.3587113866866219, 0.6: 0.3366132595640836, 0.7: 0.3308098581783994, 0.8: 0.30668543275393273, 0.9: 0.25586642990003705, 1.0: 0.19367291099489414, 1.1: 0.140509402319722, 1.2: 0.10488567232900191, 1.3: 0.0812364807836786, 1.4: 0.06121836858943995, 1.5: 0.04557161888610211, 1.6: 0.04443007843526534, 1.7: 0.06561253470317936, 1.8: 0.10221747955819366, 1.9: 0.13227159026992924, 2.0: 0.1329210655233114, 2.1: 0.09855090483145529, 2.2: 0.047674796368554234, 2.3: 0.011693423760589221, 2.4: 0.01310025279069152, 2.5: 0.0495034260404082, 2.6: 0.09560474983328712, 2.7: 0.12138942900862644, 2.8: 0.11246525862569816, 2.9: 0.07755055252846635, 3.0: 0.038910997546266854, 3.1: 0.014928113082244105, 3.2: 0.009054006524575536, 3.3: 0.012589849241910476, 3.4: 0.01647299874239502, 3.5: 0.0202763181419398, 3.6: 0.030614666162991806, 3.7: 0.0520115464150108, 3.8: 0.08083538376318122, 3.9: 0.10999796024130568, 4.0: 0.14106762557191635}

[jcmgray]

Hi @sajjan02purdue, I don't see anything obvious but sadly its too long an example for me to process what really what is going on or understand.

I would recommend trying to separate out the circuit specification from the simulation to first check if the generated gates are correct (i.e. try with the exact same list of gates in qiskit etc.). Then see if it is an issue with the MPS simulator or not.