# Error while calculating out-of-time correlators

**URL:** <https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603>\
**Category:** ITensor Julia Questions\
**Tags:** julia, mps, gate\_evolution\
**Created:** [April 8, 2026, 12:06pm UTC](https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603 "2026-04-08T12:06:52Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![asing](https://avatars.discourse-cdn.com/v4/letter/a/c68b51/32.png) [@asing](https://itensor.discourse.group/u/asing)\
**Post date:** [April 8, 2026, 12:06pm UTC](https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603/1 "2026-04-08T12:06:52Z")

</div>

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.

```auto
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.

```julia

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

```

---

<div class="post-metadata">

**Author:** ![corbett5](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/corbett5/32/886_2.png) [@corbett5](https://itensor.discourse.group/u/corbett5)\
**Post date:** [April 8, 2026, 4:25pm UTC](https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603/2 "2026-04-08T16:25:25Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![asing](https://avatars.discourse-cdn.com/v4/letter/a/c68b51/32.png) [@asing](https://itensor.discourse.group/u/asing)\
**Post date:** [April 8, 2026, 5:06pm UTC](https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603/3 "2026-04-08T17:06:28Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![corbett5](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/corbett5/32/886_2.png) [@corbett5](https://itensor.discourse.group/u/corbett5)\
**Post date:** [April 8, 2026, 10:04pm UTC](https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603/4 "2026-04-08T22:04:50Z")

</div>

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`.

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

```

Should be

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

```

with this change I get

```plaintext
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

```

---

<div class="post-metadata">

**Author:** ![Gauthameshwar](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/gauthameshwar/32/844_2.png) [@Gauthameshwar](https://itensor.discourse.group/u/Gauthameshwar)\
**Post date:** [April 10, 2026, 11:18am UTC](https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603/5 "2026-04-10T11:18:38Z")

</div>

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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![asing](https://avatars.discourse-cdn.com/v4/letter/a/c68b51/32.png) [@asing](https://itensor.discourse.group/u/asing)\
**Post date:** [April 10, 2026, 11:34am UTC](https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603/6 "2026-04-10T11:34:50Z")

</div>

Thank you!  
It was a conceptual error.

---

<div class="post-metadata">

**Author:** ![system](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/d3072b13e047cd06df7f3981547c0917d940dfa4.png) [@system](https://itensor.discourse.group/u/system)\
**Post date:** [April 20, 2026, 11:35am UTC](https://itensor.discourse.group/t/error-while-calculating-out-of-time-correlators/2603/7 "2026-04-20T11:35:11Z")

</div>

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.