Thanks for the extra details.

So actually youâ€™re in luck and this is the simplest kind of MPO to construct. You could still do it using OpSum but it would be overkill in this case. All you have to do in this case is the following. Make an MPO for your N qubits and a collection of site indices you are going to use

```
N = 10
s = siteinds("Qubit",N)
M = MPO(s)
```

then just set the MPO tensors one by one to the values you want, following the ITensor convention that the jâ€™th MPO tensor has the indices `s[j], s[j]'`

.

So like this:

```
for j=1:N
T = ITensor(s[j],s[j]')
T[1,1] = ...
T[2,1] = ...
T[1,2] = ...
T[2,2] = ...
M[j] = T
end
```

You could either work out the matrix elements of T ahead of time and then set them one by one like I am suggesting in the code snippet above. Or you could use some built in ITensor facilities for making each MPO tensor. What I mean is, say you want to make the operator

3|x+\rangle\langle x\!+\!| - I

You could do this with the code

```
T = 3*state("X+",s[j]')*state("X+",s[j]) - op("I",s[j])
```

Here I havenâ€™t included the Uâ€™s explicity, but based on what you said they are Pauli matrices so you could work out what U^\dagger_j |b_j\rangle is in advance without numerically acting U^\dagger_j onto a state.

For a list of predefined states and operators for â€śQubitâ€ť sites please see this page:

https://itensor.github.io/ITensors.jl/dev/IncludedSiteTypes.html#â€śQubitâ€ť-SiteType

You can also define your own custom overloads of the `state`

and `op`

functions if you wish:

https://itensor.github.io/ITensors.jl/dev/examples/Physics.html