create particular mixed state

ITensor Julia Questions
1

Hello everyone,

I am new to using ITensors so forgive me if I ask something simple.

I want to use purification to construct a thermal density matrix. In my specific case I think I can bypass the imaginary time evolution because the state is very close to a product state. In particular let H_{x,x+1} = H_x \otimes H_{x+1} be the local (two sites) Hilbert space. I have a spin 1/2 so H_x = \text{span}\{\ket{\uparrow}, \ket{\downarrow}\} for all x. Let us say we have 2N sites in total so that on odd sites there is the physical degree of freedom while on even sites I put the ancilla. I have to construct the following density matrix (apart an overall normalisation factor)

\rho \propto \prod_{x=1}^{N-1}\left(a\ket{\uparrow_x \uparrow_{x+1}} + b\ket{\downarrow_x \downarrow_{x+1}}\right)\left(a^*\bra{\uparrow_x \uparrow_{x+1}} + b^*\bra{\downarrow_x \downarrow_{x+1}}\right)

where a,b are two complex constants. I would expect this to be relatively easy since constructing an MPS representation of a product state is easy using the function MPS. However I am having some troubles because I am not able to construct the elementary factors

\ket{\psi_x}=a\ket{\uparrow_x \uparrow_{x+1}} + b\ket{\downarrow_x \downarrow_{x+1}}

and combine them into a density matrix without having Julia complaining about memory issues. I will really appreciate help and some explanation of it.

Eventually then, I will want to do the time evolution of this state and compute

\rho(t) = e^{i H t} \rho e^{-i H t}

and to do this I will construct appropriate ITensor gates and use the apply function. I was able to run the time evolution for simple product states as Neel states \ket{\uparrow \downarrow\dots \uparrow\downarrow} so the rest of my code works. The only problem I have is to construct the initial state above.

I really thank you for the help!

Giuseppe

2

Hello Giuseppe,

Can you please post a snippet of the Julia code you use to construct this initial state. If possible, can you please also show what is thrown by Julia regarding the memory issues you describe.

Thanks,
Karl

3

Hello Karl thanks for your reply. In the end I have solved my problem. Specifically, I had to create the following “thermal” density matrix

\rho = \rho_{\mu_L}\otimes \rho_{\mu_R}

where

\rho_\mu = \frac{1}{Z}\prod_{x=1}^{N-1}(e^{-\mu/2} \ket{\uparrow_{x} \uparrow_{x+1}} + e^{\mu/2}\ket{\downarrow_x\downarrow_{x+1}})(e^{-\mu/2} \bra{\uparrow_{x} \uparrow_{x+1}} + e^{\mu/2}\bra{\downarrow_x\downarrow_{x+1}})

with Z= 2 \cosh(\mu/2) . In particular choosing N=n/2 with n even, we are looking for a state with different (average) magnetisation on the left and the right of the middle of the chain. This state is also known as domain wall, or variations of this name. My approach has been the following:

function initial_state(Nphys::Int64, μL::Float64, μR::Float64)
      N = 2Nphys
      sites = siteinds("S=1/2", N)
      state = MPS(sites)

      #create the tensors on the two sites in the first half
      for x in 1:2:Nphys-1
            # copy the tensors
            t1 = state[x]
            t2  = state[x+1]
            # create new indices

            x_physical1 = sites[x] # first physical index
            i_virtual = Index(2, tags(commonind(t1, t2))) # bond connecting the two tensors 
            x_physical2 = sites[x+1] #second physical index (the ancilla)

            # construct the tensors

            A = ITensor(x_physical1, i_virtual)
            B = ITensor(x_physical2, i_virtual)
            # x_physical = 1 is up, x_physical = 0 is down
            A[x_physical1 => 1, i_virtual => 1] = 1/(2cosh(μL)) * exp(-μL) |> sqrt
            A[x_physical1 => 1, i_virtual => 2] = 1/(2cosh(μL)) * exp(-μL) |> sqrt
            A[x_physical1 => 2, i_virtual => 1] = 1/(2cosh(μL)) * exp(μL) |> sqrt
            A[x_physical1 => 2, i_virtual => 2] = 1/(2cosh(μL)) * exp(μL) |> sqrt

            B[x_physical2 => 1, i_virtual => 1] = 1/(2cosh(μL)) * exp(-μL) |> sqrt
            B[x_physical2 => 1, i_virtual => 2] = 1/(2cosh(μL)) * exp(-μL) |> sqrt
            B[x_physical2 => 2, i_virtual => 1] = 1/(2cosh(μL)) * exp(μL) |> sqrt
            B[x_physical2 => 2, i_virtual => 2] = 1/(2cosh(μL)) * exp(μL) |> sqrt

            # put these tensor into the MPS

            state[x], state[x+1] = A, B
            sites[x], sites[x+1] = x_physical1, x_physical2
      end

      print("Finished first half...\n")
      # second half
      for x in Nphys+1:2:N-1
            # copy the tensors
            t1 = state[x]
            t2  = state[x+1]
            # create new indices

            x_physical1 = sites[x] # first physical index
            i_virtual = Index(2, tags(commonind(t1, t2))) # bond connecting the two tensors 
            x_physical2 = sites[x+1] #second physical index (the ancilla)

            # construct the tensors

            A = ITensor(x_physical1, i_virtual)
            B = ITensor(x_physical2, i_virtual)
            # x_physical = 1 is up, x_physical = 0 is down
            A[x_physical1 => 1, i_virtual => 1] = 1/(2cosh(μR)) * exp(-μR) |> sqrt
            A[x_physical1 => 1, i_virtual => 2] = 1/(2cosh(μR)) * exp(-μR) |> sqrt
            A[x_physical1 => 2, i_virtual => 1] = 1/(2cosh(μR)) * exp(μR) |> sqrt
            A[x_physical1 => 2, i_virtual => 2] = 1/(2cosh(μR)) * exp(μR) |> sqrt

            B[x_physical2 => 1, i_virtual => 1] = 1/(2cosh(μR)) * exp(-μR) |> sqrt
            B[x_physical2 => 1, i_virtual => 2] = 1/(2cosh(μR)) * exp(-μR) |> sqrt
            B[x_physical2 => 2, i_virtual => 1] = 1/(2cosh(μR)) * exp(μR) |> sqrt
            B[x_physical2 => 2, i_virtual => 2] = 1/(2cosh(μR)) * exp(μR) |> sqrt

            # put these tensor into the MPS

            state[x], state[x+1] = A, B
            sites[x], sites[x+1] = x_physical1, x_physical2
      end
      print("Finished second half...\n")
      sites, state

end

You can check that computing the expectation value of S^z_x results in the correct thing. Maybe there’s a smarter way though but at the moment the function looks very efficient…

Thank you!

1 Like
4

Giuseppe,

Great! I am so glad that you were able to figure this out and that the error was unrelated to the ITensors software. I took a minute to help write a version of this initial_state function which reduces arithmetic which should help when you have many sites, you can find it below. I verified that the numerical state out is equal to the state produced by your code.

All the best,
Karl

function create_states(state::MPS, sites, matrix, first, last)
  for x in first:2:last
    # copy the tensors
    t1 = state[x]
    t2  = state[x+1]
    # create new indices

    x_physical1 = sites[x] # first physical index
    i_virtual = Index(2, tags(commonind(t1, t2))) # bond connecting the two tensors 
    x_physical2 = sites[x+1] #second physical index (the ancilla)

    # construct the tensors

    A = ITensor(matrix, (x_physical1, i_virtual))
    B = ITensor(matrix, (x_physical2, i_virtual))

    # put these tensor into the MPS

    state[x], state[x+1] = A, B
    sites[x], sites[x+1] = x_physical1, x_physical2
  end
end

function initial_state(Nphys::Int64, μL::Float64, μR::Float64)
  N = 2Nphys
  sites = siteinds("S=1/2", N)
  state = MPS(sites)

  #create the tensors on the two sites in the first half
  cosuL = 1/(2cosh(μL))
  expmuL = exp(-μL)
  expuL = exp(μL)

  val1 = cosuL * expmuL |> sqrt
  val2 = cosuL * expuL |> sqrt

  matrix = [val1, val2, val1, val2]
  create_states(state, sites, matrix, 1, Nphys -1)

  print("Finished first half...\n")
  cosuR = 1/(2cosh(μR))
  expmuR = exp(-μR)
  expuR = exp(μR)
  val1 = cosuR * expmuR |> sqrt
  val2 = cosuR * expuR |> sqrt

  matrix = [val1, val2, val1, val2]
  create_states(state, sites, matrix, Nphys+1, N-1)

  print("Finished second half...\n")
  sites, state

end
2 Likes