Complex values of time dependent single site spin correlations

ITensor Julia Questions
,
1

Hello!

I am trying to calculate \langle S^{x}_1(t)S^{x}_1(0)\rangle for a randomMPS initial state |\psi \rangle but the output of the correlations have imaginary parts of the order 10^{-3}.

I am using the Julia version of the code but have more experience with the C++ version and so its possible that I do not understand all the functions I am using correctly. My current approach is to obtain |\psi(t)\rangle = U(t) |\psi\rangle and | \psi_x(t) \rangle= U(t) S_x |\psi\rangle and then calculate \langle \psi(t) | S_x | \psi_x (t)\rangle which should equal the correlator I am seeking.

I am using TEBD for the time evolution generated by two different Hamiltonians. I obtain the time evolved states by incrementing by some small dt until t = nT where I measure the overlap. To get the two states I do:

    psi = randomMPS(sites,bond_dim)
    psi_x = deepcopy(psi)
    
    #apply Sx to psi_x 
    opp = op("Sx",sites[1])
    A_1_x = opp * psi_x[1]
    noprime!(A_1_x)
    psi_x[1] = A_1_x

I initially check the overlap at t=0 to ensure things are fine

val = inner(psi',Sx,psi_x)
println("<Sx_1(0)Sx_1(0)> = $val")

It gives the correct value of 0.25. I then time evolve using the second order trotter decomposition

U(\delta) = e^{-i H_{even} \delta/2} e^{-i H_{odd} \delta} e^{-i H_{even} \delta/2}

and normalize after each step (unsure why I do this but it was included in the TEBD code example). Once the evolutionis done at an n where t = nT I then calculate the overlap between the states with Sx as an MPO.

#perform time evolution
    for n in 1:cycles
        #first apply H1 with bricklaying pattern for T/2
        c = 1
        for t in 0.0:tau:T/2
            #println("$((n-1)*T + t)")
            psi_x = apply(gates1_e,psi_x; cutoff)
            psi_x = apply(gates1_o,psi_x; cutoff)
            psi_x = apply(gates1_e,psi_x; cutoff)
            normalize!(psi_x)
            psi = apply(gates1_e,psi; cutoff)
            psi = apply(gates1_o,psi; cutoff)
            psi = apply(gates1_e,psi; cutoff)
            normalize!(psi)
        end
        #Apply H2 afterwords for T/2
        for t in 0.0:tau:T/2
            psi_x = apply(gates2,psi_x; cutoff)
            normalize!(psi_x)
            psi = apply(gates2,psi; cutoff)
            normalize!(psi)
        end
        #get overlap with initial state:  
        val = inner(psi',Sx,psi_x)/2
        push!(overlaps,real(val))
        println("<Sx_1($(n)T) * Sx_1(0)> = $val")

    end

The results for N = 8 and bond dim=60 (for speed of debugging) give:

<Sx_1(0)Sx_1(0)> = 0.25000000000000006
<Sx_1(1T) * Sx_1(0)> = 0.21372936110782428 - 0.02270611893409933im
<Sx_1(2T) * Sx_1(0)> = 0.12311194059048045 - 0.06449038038510997im
<Sx_1(3T) * Sx_1(0)> = 0.015579761828750855 - 0.10712389554900748im
<Sx_1(4T) * Sx_1(0)> = -0.07714945391656257 - 0.11487553717297402im
<Sx_1(5T) * Sx_1(0)> = -0.13803104880264594 - 0.07710401302521291im
<Sx_1(6T) * Sx_1(0)> = -0.15916753080751667 - 0.010361880398627841im
<Sx_1(7T) * Sx_1(0)> = -0.13946923498825858 + 0.05843952708266153im
<Sx_1(8T) * Sx_1(0)> = -0.08687407222840543 + 0.10668196023003335im
<Sx_1(9T) * Sx_1(0)> = -0.01785099025957232 + 0.12293533638047521im
<Sx_1(10T) * Sx_1(0)> = 0.047433252646377363 + 0.10736005379295518im

I tried other ways of getting the inner product by doing it by hand with the individual matrices but it yields similar results. The values will be compared with ED results later on, I just dont know why im getting imaginary parts to these correlations as they should be real (unless I am mistaken).

Any help would be great!
Sebastien

2

It’s a good question, but I suppose since the real-time evolution involves imaginary numbers, that will of course make the states complex-valued, so then in general observables can have imaginary parts. Since your inner product involves two (technically) different states, then there’s no symmetry of the inner product that would force the imaginary part to be zero. (I’m just thinking through whether this is a bug or a possible thing – I think it’s totally possible and mathematically allowed.)

So then the next things to check would be: what happens if you do your calculations to higher accuracy i.e. use a somewhat smaller cutoff or a smaller time step tau? Do the imaginary parts go down?

3

Hi Miles, thanks for the quick reply.

I agree that it looks mathematically possible to have complex values, for some reason Im slightly surprised as its just a spin-spin correlation function.

I profiled several cases with different parameters. In general, it would seem that slowly increasing \chi (from 50, 70, 80,100, 200) while decreasing the cutoff (from 10^{-8} to 10^{-10}) and the time step \tau (from 0.1, 0.01,0.005, 0.001) has no effect on the magnitude of the imaginary part. The imaginary part vary more or less between orders of 10^{-2} \text{ and } 10^{-3}, while the real parts remain sensical but are consistently larger than the imaginary part - typically always at scales of 10^{-1} . I’ve never used time evolution evolution with MPS’s and MPO’s so I dont know if this small imaginary part is to due to the errors accrued at each step and I am trying to study too far in time or that I did something else wrong. I will include the last printout of results for N = 8, \chi = 300, \tau = 0.001 and cutoff = 10^{-10} for 10 cycles. Some imaginary part are smaller than what I claimed above but others are of the same order.

<Sx_1(0)Sx_1(0)> = 0.2500000000000001
<Sx_1(1T) * Sx_1(0)> = 0.2096352776653267 - 0.002072090428430637im
<Sx_1(2T) * Sx_1(0)> = 0.11886952980874779 - 0.006964053883293589im
<Sx_1(3T) * Sx_1(0)> = 0.019216187442766257 - 0.01145594845359801im
<Sx_1(4T) * Sx_1(0)> = -0.06405043721444077 - 0.013252421054170353im
<Sx_1(5T) * Sx_1(0)> = -0.12035980414857056 - 0.012438981128730134im
<Sx_1(6T) * Sx_1(0)> = -0.14922265410864494 - 0.010181901669937217im
<Sx_1(7T) * Sx_1(0)> = -0.15214553105493317 - 0.0080472314030757im
<Sx_1(8T) * Sx_1(0)> = -0.1302834235849751 - 0.006516715074302525im
<Sx_1(9T) * Sx_1(0)> = -0.08639714354327803 - 0.00424560130572748im
<Sx_1(10T) * Sx_1(0)> = -0.02710903403262373 + 0.000982732658052376im

Are there any good measures for determining the accuracy of a time evolved state like this? I still am unsure about good metrics to track that can inform me about the quality of the evolution in time. Some preliminary ED results also might have an imaginary part which I will also be checking.

Best,
Sebastien

4

Re-ran once more for \chi = 400 with the same other parameters posted above in the last printout:

<Sx_1(1T) * Sx_1(0)> = 0.21490454476528975 - 0.025244828021217452im
<Sx_1(2T) * Sx_1(0)> = 0.12297362104960749 - 0.055278756913562294im
<Sx_1(3T) * Sx_1(0)> = 0.006155582644072806 - 0.06991997303735265im
<Sx_1(4T) * Sx_1(0)> = -0.09090175265124087 - 0.049363880809939135im
<Sx_1(5T) * Sx_1(0)> = -0.1391883654795017 - 0.0007733187375506857im
<Sx_1(6T) * Sx_1(0)> = -0.1358487308698457 + 0.04982855636268106im
<Sx_1(7T) * Sx_1(0)> = -0.0944724484073148 + 0.07887102052358512im
<Sx_1(8T) * Sx_1(0)> = -0.03522666965367462 + 0.07687277760441888im
<Sx_1(9T) * Sx_1(0)> = 0.022292626501502988 + 0.04876632014498093im
<Sx_1(10T) * Sx_1(0)> = 0.06370853663598953 + 0.008047497085259928im

I do need to use a randomMPS initial state for this but it seems like it can have very mixed results as a consequence.

5

It might just be that these imaginary parts are really there. Like e.g. just as a pedagogical example, if the initial state was the ground state of your Hamiltonian, then the time evolution operator acting in the left part of your correlator (between the left state and left Sx, the one at time t) just acts on your initial state to give an overall time-dependent phase of e^{iE_0 t} (I hope I got the sign correct there). Whereas the other time-evolution operator is in between two Sx operator so does something more complicated and doesn’t give an opposite phase. So there could be some general imaginary part.

About testing, ED is a great way. Another way could be to time evolve your MPS by some amount forward then time evolve it next by the same amount backward, back to time zero. Now overlap with your initial state (save a copy of your initial state to do this) and see how close the fidelity is to 1.0. This could let you determine what might be some good values of cutoff, maxdim, and time step to take for the accuracy you want.

6

Sounds good,

I will take a look into the physics and see if If I am missing a factor somewhere. Thanks for the info on the fidelity, that’s a good tip!

Thank you,
Sebastien