Error while calculating out-of-time correlators

ITensor Julia Questions
, ,
1

Hi!

In the process of calculating OTOC’s for a spin chain mode, I came across a strange error when I perform backward time evolution.
In the following code, I am performing

C_{i1}(i,t) = \langle \psi_{init} | S^z_i(t) S^z(0) |\psi_{init}\rangle

I am doing it in two different ways

  1. \psi_R (t) = S^z_i e^{-iHt} S^z_0|\psi_{init \rangle}
    \psi_L (t) = e^{-iHt} |\psi_{init \rangle}
    C[i,t] = \langle \psi_L(t) , \psi_R(t) \rangle
  2. \phi (t) = e^{+iHt }S^z_i e^{-iHt} S^z_0|\psi_{init} \rangle
    F[i,t] = \langle \psi_{init} , \phi(t) \rangle

Logically, both the methods should give the same result but I notice discrepancy in the result.

Starting TEBD time evolution...
Steps: 40  |  dt = 0.05  |  t_end = 2.0
------------------------------------------------------------
t =  0.50  |  C(L/2,t) =  0.6849 + -0.0000i  |  F(L/2)=-0.4024+-0.4702i | abs(F) =  0.6189 | maxdim = 4
t =  1.00  |  C(L/2,t) =  0.5726 + -0.0000i  |  F(L/2)=-0.5041+ 0.1552i | abs(F) =  0.5274 | maxdim = 6
t =  1.50  |  C(L/2,t) =  0.5867 +  0.0000i  |  F(L/2)= 0.0110+ 0.5736i | abs(F) =  0.5737 | maxdim = 8
t =  2.00  |  C(L/2,t) =  0.6443 +  0.0000i  |  F(L/2)= 0.7891+ 0.0646i | abs(F) =  0.7918 | maxdim = 8

I have checked that the states are normalized at each time step and after application of certain operators.


using ITensors
using ITensorMPS
using Printf
using DelimitedFiles
using Plots

# ── parameters ────────────────────────────────────────────────
const L      = 6       # chain length
const J      = 1.0      # Ising coupling (ZZ)
const hx     = 1.05    # transverse field (x)
const hz     = 0.5      # longitudinal field (z)
const j_ref  = 1        # reference site for correlator (1-indexed)

const dt     = 0.05     # time step
const t_end  = 2.0     # total evolution time
const maxdim = 200      # max MPS bond dimension
const cutoff = 1e-10    # SVD truncation cutoff

function make_gates(sites, J, hz, hx, dt)
    L     = length(sites)
    gates = ITensor[]
   
    for j in 1:L-1
        s1, s2 = sites[j], sites[j+1]
        hj = ITensor()
        #hj = emptyITensor(s1, s2)
        # Construct the two-site gate for the interaction term
        if j == 1
            # For the first and last bond, we only have one site with the interaction term
             hj =  -J*4.0* op("Sz", s1) * op("Sz", s2) 
             hj += -hx*2.0* op("Sx", s1) * op("Id", s2) 
             hj += -hz*2.0* op("Sz", s1) * op("Id", s2)
             hj += -hx/2*2.0* op("Id", s1) * op("Sx", s2)
             hj += -hz/2*2.0* op("Id", s1) * op("Sz", s2)
        elseif j == L-1
             hj =  -J*4.0* op("Sz", s1) * op("Sz", s2) 
             hj += -hx/2*2.0* op("Sx", s1) * op("Id", s2)
             hj += -hz/2*2.0* op("Sz", s1) * op("Id", s2)
             hj += -hx*2.0* op("Id", s1) * op("Sx", s2)
             hj += -hz*2.0* op("Id", s1) * op("Sz", s2)
        else
            # For the middle bonds, we have both sites with the interaction term
             hj =  -J*4.0* op("Sz", s1) * op("Sz", s2) 
             hj += -hx/2*2.0* op("Sx", s1) * op("Id", s2) 
             hj += -hz/2*2.0* op("Sz", s1) * op("Id", s2)
             hj += -hx/2*2.0* op("Id", s1) * op("Sx", s2)
             hj += -hz/2*2.0* op("Id", s1) * op("Sz", s2)
            
        end
        # gate: e^{-i hj dt/2}
        Gj = exp(-im * dt/2 * hj)
        push!(gates, Gj)
    end
     append!(gates, reverse(gates))  # Add the gates in reverse order for the second half of the step
    return gates   # length (2*(L-1))
end

let 
n_steps = Int(round(t_end / dt))
times   = collect(0.0:dt:t_end)

# ── sites ─────────────────────────────────────────────────────
sites = siteinds("S=1/2", L; conserve_qns=false)
gates = make_gates(sites, J, hz, hx, dt)
gates_dag = make_gates(sites, J, hz, hx, -dt)  # gates for backward evolution (dagger)
println("Gates built: $(length(gates)) two-site gates")
println("Gates built: $(length(gates_dag)) two-site gates")


#── TEBD time evolution ───────────────────────────────────────

psi_init = MPS(sites, "Up") # Start from a simple product state (all spins down)

# ── apply Sz at reference site to initial state ────────────────
Sz_ref    = 2.0*op("Sz", sites, j_ref)
psi_R     = apply(Sz_ref, psi_init)     # |phi_R> = Sz_{j_ref}|psi_init>

C_it  = zeros(ComplexF64, L, length(times)) 
F_it  = zeros(ComplexF64, L, length(times))   # OTOC, <psi| Sz_i(t) Sz_0 Sz_i(t) Sz_0 |psi>

# ── measure at t=0 ────────────────────────────────────────────

for i in 1:L
    Sz_i    = 2.0*op("Sz", sites, i)

    Sz_i_psi_R   = apply(Sz_i, psi_R)
    C_it[i, 1]   = inner(psi_init, Sz_i_psi_R)

    F_it[i, 1]   = C_it[i, 1] 
end

# ── time evolution loop ───────────────────────────────────────
println("Starting TEBD time evolution...")
println("Steps: $n_steps  |  dt = $dt  |  t_end = $t_end")
println("-"^60)

psi_L = copy(psi_init)   # |psi(t)>

for step in 1:n_steps
    t = step * dt
    # evolve both states by one dt
    psi_L = apply(gates, psi_L; maxdim, cutoff) # e^{-iHt} |psi(t)>
    psi_R = apply(gates, psi_R; maxdim, cutoff) # e^{-iHt} [Sz_{j_ref} |psi(t)>]
    for i in 1:L
        Sz_i = 2.0*op("Sz", sites, i)

        Sz_i_psi_R = apply(Sz_i, psi_R)
        C_it[i, step+1] = inner(psi_L, Sz_i_psi_R) # <psi(t)| Sz_i [e^{-iHt} Sz_{j_ref} |psi(t)>]

        phi = apply(Sz_i, psi_R)        # Sz_i [e^{-iHt} Sz_0 |psi>]
        phi = apply(gates_dag, phi; maxdim, cutoff) 

        F_it[i, step+1] = inner(psi_init, phi)  # <psi_init| e^{+iHt} Sz_i [e^{-iHt} Sz_0 |psi_init>]
        #F_it[i, step+1] = inner(psi_L, phi)  # <psi(t)| Sz_i [e^{-iHt} Sz_0 |psi>]
    end

     if step % 10 == 0
        @printf("t = %5.2f  |  C(L/2,t) = %7.4f + %7.4fi  |  F(L/2)=%7.4f+%7.4fi | abs(F) = %7.4f | maxdim = %d\n",
                t, 
                real(C_it[L÷2, step+1]), imag(C_it[L÷2, step+1]),
                real(F_it[L÷2, step+1]), imag(F_it[L÷2, step+1]),
                abs(F_it[L÷2, step+1]),
                maxlinkdim(psi_L))
    end
end
end
2

It’s more accurate to evolve two states a time t (method 1) than it is to evolve a single state a time 2t.

If this is the case I would expect the two methods to agree at earlier time steps. You could also try increasing the cutoff, since your MPS bond dimension remains very small at all times.

3

Hi corbett5,
Thank you for your message. It does make sense to follow method 1, since the second method will accumulate more errors. However, as mentioned in the post I found this error while performing the calculation for OTOC where I need to perform minimum two time evolutions (forward and backward).

I have tried increasing the cutoff but the answer has not changed since the system size is too small.

4

Ok I think I know why the error is occurring. You are only backwards evolving by dt, but you would really need to backwards evolve another step * dt.

phi = apply(gates_dag, phi; maxdim, cutoff)

Should be

for _ in 1:step
  phi = apply(gates_dag, phi; maxdim, cutoff) 
end

with this change I get

Steps: 40  |  dt = 0.05  |  t_end = 2.0
------------------------------------------------------------
t =  0.50  |  C(L/2,t) =  0.6849 +  0.0000i  |  F(L/2)= 0.6849+-0.0000i | abs(F) =  0.6849 | maxdim = 8
t =  1.00  |  C(L/2,t) =  0.5726 +  0.0000i  |  F(L/2)= 0.5726+-0.0000i | abs(F) =  0.5726 | maxdim = 8
t =  1.50  |  C(L/2,t) =  0.5867 + -0.0000i  |  F(L/2)= 0.5867+-0.0000i | abs(F) =  0.5867 | maxdim = 8
t =  2.00  |  C(L/2,t) =  0.6443 +  0.0000i  |  F(L/2)= 0.6443+-0.0000i | abs(F) =  0.6443 | maxdim = 8
1 Like
5

That makes total sense. I tried it as well, and it worked. To avoid bugs like these, I usually define a helper function beforehand and call it within the loops so that these elusive “nested loops” don’t cause unnecessary bugs. For example, this would make the backward evolution obvious in each loop:

function evolve_state(psi, gates, steps; maxdim, cutoff)
    for _ in 1:steps
        psi = apply(gates, psi; maxdim=maxdim, cutoff=cutoff)
    end
    return psi
end

for step in 1:n_steps 
    psi_L = apply(gates, psi_L; maxdim, cutoff) 
    psi_R = apply(gates, psi_R; maxdim, cutoff) 
    
    for i in 1:L 
        phi = apply(op("Sz", sites, i), psi_R) 
        # Explicitly evolve 'phi' all the way back to the start 
        phi_at_zero = evolve_state(phi, gates_dag, step; maxdim, cutoff) 
  
        F_it[i, step+1] = inner(psi_init, phi_at_zero) 
    end 
end
6

Thank you!
It was a conceptual error.