General Question about MPO's and the apply function

Hey all!
I have some urgent questions about ITensor.
The process I am trying to calculate is built from four main part:

  1. Building a density matrix.
  2. Building a long range interaction gates.
  3. Develop the density matrix in unitary time using the “apply” (TEBD) .
  4. Calculate Tr[\rho S^+], Where is: S^+=I \otimes ... \otimes I \otimes S^+ .

And unfortunately I have questions regarding to some of them.

About the first part:

The way I thought to construct the density matrix is using matrices that represents the states on each site, and the making it an ITensor object of the form of MPO.
Example:

state=1/2 *[1 1 ; 1 1]
rho= state
for j=2:n
   rho=kron(rho, state)
end
rho_t=ITensor(rho, site,site')
r=MPO(rho_t, site)

the problem is that as it seems this “r” will have bond dimensions greater than 1, how is it possible ? because I am a little bit confused about the bond dimension because I thought it should represent the amount of entanglement in the bipartite of site and there is non.
Also, why their is “Link,n=1” and sometimes “Link,l=1” ?

MPO
[1] ((dim=2|id=641|"Qubit,Site,n=1"), (dim=2|id=641|"Qubit,Site,n=1")', (dim=4|id=32|"Link,n=1"))
[2] ((dim=4|id=32|"Link,n=1"), (dim=2|id=251|"Qubit,Site,n=2"), (dim=2|id=251|"Qubit,Site,n=2")', (dim=4|id=521|"Link,n=2"))
[3] ((dim=4|id=521|"Link,n=2"), (dim=2|id=161|"Qubit,Site,n=3"), (dim=2|id=161|"Qubit,Site,n=3")')

This is the output in the case of n=3.

about the second part:
Is there any way to use long range interactions to simulate TEBD ? I know that you can use sweeps, but is it accurate ?

Now I will ask about the fourth part.

for the example I wrote above Tr[\rho S^+]= 0.5. The way that I am calculating this is using the apply:

tr(apply(r,S_plus))

Does this is the correct way on doing so? Because it seems that something is going wrong.
For example if I apply only interactions from site 1:n-1 and then calculating the trace I am getting the the trace is decreasing over time (and this is not possible) and it happens even non- interacting systems (if needed I can provide a simple example).
The way I am developing in time is using the “apply”.

Example for gates:

    sz=[1 0; 0 -1]
    gates=ITensor[]
    for i=1:n-1
        a= sz
        hin = ITensor(a, (site[i]), (site[i])')
        Hin = exp((tau / 2) * hin)
        push!(gates, Hin)
    end
    append!(gates, reverse(gates))

Is there any reason why entanglement is building in a system like that ?

I hope that someone can help me with those issues.

Can anyone help me with that issue please ?

I also tried to check if the density matrix I made is equivalent to that one:

dens=MPO(site , n->state)

And the two objects aren’t equivalent.