Hi all, I’ve seen a couple of questions regarding circuits with mixed state initialization, but I have some further questions. I would like to use `QubitDensityMatrix`

to implement my mixed state as input however, I realised a couple of problems that I could not solve. It seems like `QubitDensityMatrix`

is only available for “default.mixed” devices. I’ve seen here that mixed states are also available in “cirq.mixedsimulator” device but `qml.device("cirq.mixedsimulator", wires=2).operations`

does not have `QubitDensityMatrix`

. I want to use it in TensorFlow because, for some reason, I can not take the gradient of my expectation value if I use `scipy.linalg.expm`

to construct my hamiltonian. Here is a simplified example;

Imagine that I have a trainable function which creates my Hamiltonian, and I need to take the gradient of the expectation value wrt the parameters of this function and quantum network. Here is some code

```
from scipy.linalg import expm
dev = qml.device("default.mixed", wires = 2)
rho = np.zeros((2**2,2**2), dtype=np.complex128)
rho[0,0] = 1.
state = np.random.uniform(0.1, 1,(2**2, 1))
state /= np.linalg.norm(state)
energy = -3.
phi = np.random.uniform(0, 1,(2,))
@qml.qnode(dev)
def circuit(rho: np.ndarray, phi: np.ndarray):
qml.QubitDensityMatrix(rho, wires=range(2))
qml.RY(phi[0], 0)
qml.RY(phi[1], 1)
qml.CNOT(wires=[0,1])
return qml.density_matrix(wires=range(2))
def cost(rho, phi,state, energy):
Hamiltonian = energy * (state @ np.conj(state).T)
density = circuit(rho, phi)
return np.real(np.trace(density @ expm(-Hamiltonian)))
grads = qml.grad(cost)
g=grads(rho, phi, state, energy)
```

Here since my Hamiltonian is an ArrayBox, this code crushes. But without ` expm(-Hamiltonian)`

, I get a result.

The reason why I didn’t use `qml.expval(qml.Hermitian(Hamiltonian, wires=range(2)))`

is because I again got an Arraybox-related issue where this time, somewhere in the code, it’s trying to take conjugate of the matrix which does not exist since its an ArrayBox. Note that this works if I run `cost(rho, phi, state, energy)`

without gradient. My only issue is with the gradient. So I thought the simplest workaround would be to externalize the gradients to TensorFlow.

Is there any workaround that I can use? In the worst-case scenario, is there a generic algorithm to build a gate structure that creates a mixed state input?

Thanks