How to derive thermal state by imaginary time evolution (purification)

hi
I want to calculate the correlation function \left\langle\psi\right|e^{iHt}Ae^{-iHt}B\left|\psi\right\rangle by TEBD.
To derive the thermal states, I do the imaginary time evolution on enlarged Hilbert space \left | n \right \rangle \otimes \left | \tilde{n} \right \rangle, so I apply the gates to the identity operator. The sites of ancilla space are not evolved, so I add “prime” to the ancilla sites.
After imaginary time evolution, I need to calculate the correlation at moment t during real time evolution, but the function inner(MPO, MPO, MPO) does not work, so I apply the code "inner(MPO, apply(MPO,MPO)).
But the code does not work. What is the reason?
I just transferred from C++to Julia, so the code I wrote is still quite messy. Sorry.

using ITensors
using Printf

g = 0.0075 
U = 0.1
t = 0.03 
E = -0.01
x = -5
N =10

function ITensors.op(::OpName"expτSS", ::SiteType"Electron", s1::Index, s2::Index; tau)
    if s1 == 1
        factor = 1
    else
        factor1 = 0.5
    end

    s2 == N ? factor2 = 1 : factor2 = 0.5

    h = factor1 * E * op("Ntot", s1) + factor1 * U * op("Nupdn", s1)
      + factor1 * E * op("Ntot", s2) + factor1 * U * op("Nupdn", s2) 
      + t * op("Cdagup", s1) * op("Cup", s2) + t * op("Cdagdn", s1) * op("Cdn", s2)
    return exp(tau * h)
end

function main(; N=10, cutoff=1E-8, tstep=0.1, beta_max=4.0)
    println("begin")
    onsite_E = E2

    # Make an array of 'site' indices
    s = siteinds("Electron", N; conserve_qns=false)

    # Make gates (1,2),(2,3),(3,4),...
    gates_imag = ops([("expτSS", (n, n + 1), (tau = -tstep / 2,)) for n in 1:(N - 1)], s)
    # Include gates in reverse order to complete Trotter formula
    append!(gates_imag, reverse(gates_imag))

    gates_real = ops([("expτSS", (n, n + 1), (tau = -tstep / 2,)) for n in 1:(N - 1)], s)
    append!(gates_real, reverse(gates_real))

    # Hamiltonian of Hubbard Model
    ampoHub = OpSum()
    for j in 1:(N - 1)
        j == 1 ? onsite_E = E1 : onsite_E = E2
        ampoHub += E, "Ntot", j
        ampoHub += U, "Nupdn", j
        ampoHub += "Cdagup", j, "Cup", j + 1
        ampoHub += "Cdagdn", j, "Cdn", j + 1
    end
    ampoHub += E, "Ntot", N
    ampoHub += U, "Nupdn", N
    Hub = MPO(ampoHub, s)

    # Dx Hamiltonian of Hubbard Model
    ampoDxHub = OpSum()
    ampoDxHub += √2 * g, "Ntot", 1
    DxHub = MPO(ampoDxHub, s)

    # Initial state is infinite-temperature mixed state
    #rho = MPO(s, "Id") ./ √2  

    X = MPO(s, "Id")
    Y = MPO(prime(s), "Id")
    x = convert(MPS, X)
    y = convert(MPS, Y)
    rho = outer(x, y) ./ √2  


  # Do the imaginary time evolution by applying the gates for Nsteps steps
    for beta in 0 : tstep : beta_max
        rho = apply(gates_imag, rho; cutoff)
        rho = rho / tr(rho)
    end

    orho = apply(DxHub, rho; cutoff)

    for i in 0 : 0.1 : 100
        rho = apply(gates_real, rho; cutoff)
        orho = apply(gates_real, orho; cutoff)
        correlation = inner(rho, apply(DxHub, orho))
        @printf("real = %.8f  imag = %.8f", real(correlation), imag(correlation))
    end

    return nothing

end


main()

Are you trying to evaluate \left\langle\psi\right|e^{iHt}Ae^{-iHt}B\left|\psi\right\rangle? Is |\psi\rangle a state, i.e. an MPS? And what are A and B, are you assuming you have them in the form of an MPO? I think you should be a bit clearer about what you are trying to calculate, otherwise it is difficult to help.

Thanks for your reply.
I am trying to evaluate Tr$\left ( A\left ( t \right ) B\rho \right ) , \left | \psi \right \rangle $ is a MPS.
Tr\left ( A\left ( t \right ) B\rho \right )=\frac{1}{Z}\sum_{n}\left \langle n \right |e^{iHt}A e^{-iHt}Be^{-\beta H}\left | n \right \rangle
For evaluate the correlation of finite temperature, I use the purification method, so \left | \psi \right \rangle \equiv \frac{1}{\sqrt{Z} } \sum_{n}e^{\frac{-\beta E_{n} }{2}}\left | n \right \rangle \otimes\left | \tilde{n} \right \rangle . In practical calculations, I want to use the identity operator as the initial state and then perform virtual time evolution on it.
A and B are both in the form of MPO. A=B=\hat{n}_{1\uparrow }+\hat{n}_{1\downarrow }

Something that might help you here, evan, is that it is common when implementing the purification method to represent the state \ket{\psi} as an MPS. (Actually that is where the name purification comes from as you may know, namely that one maps an operator to a state.) Then on a technical level it can be easier to implement because one is just evolving a state (MPS) rather than an MPO. So you could in that case compute correlations in the usual way that you would for an MPS.

About putting primes on the ancilla sites, if you use the MPS approach I’d recommend instead that you make the ancilla sites to be actually different sites (different indices) without primes.

1 Like