I have experimented with MPOs representing a quantum circuit.
Here’s a simple example case:
using ITensors, ITensorMPS
let
s = siteinds("Qubit",2)
mpo = MPO(s, "I")
mpo = apply(op("Rx", s[1], θ=π/5), mpo)
mpo = apply(op("Ry", s[2], θ=π/3), mpo)
psi = MPS(s, "0")
psi = apply(mpo, psi)
end
My main question concerns the exact nature of the function apply: Does the apply function used with an operator acting on an MPO cause the operator tensor to be contracted with the MPO, or what is actually happening when an operator is applied to an MPO? Such a line in the above code is for example apply(op("Rx", s[1], θ=π/5), mpo)
For an MPO M and operator O_i on site i, apply(Oᵢ, M) performs the operation O_iM, i.e. it left-multiplies the operator onto the site (in practice, it finds the site by searching the MPO for the site that matches the indices of the operator). For a single-site operator it just performs a contraction of O_i with the M_i, the MPO tensor at site i, and also does some prime level manipulation to that O_iM has the same prime levels and indices as M. For a multi-site operator, it first applies swap gates to make the sites adjacent, applies the gates over the sites exactly, and then performs a truncation of the MPO according to the specified truncation parameters.
If you want to apply the operator and its conjugate, i.e. perform the operation O_iMO_i^{\dagger}, you can use apply(Oᵢ, M; apply_dag=true).
I see that the docstring for product/apply is a little bit vague on these points, we should make it more specific.