Issue with QFIM for time evolved mixed state

Thanks @CatalinaAlbornoz . There is a smaller version (no underlying physical meaning). For me it seems the applied kind of diff_method causes the error for trotter steps being set greater than 1.

dev = qml.device("default.mixed", wires=3)

# Parameter values
theta_values = pnp.array([np.pi / 4], requires_grad=True)
TROTTER_STEPS = 2    # number of Trotter steps for the approximation

@qml.qnode(dev, diff_method="parameter-shift")
def circuit(theta_values):
    
    # Construct Hamiltonians 
    H_Z = qml.Hamiltonian([theta_values[0]], [qml.PauliZ(0)])
        
    # Apply time evolution: exp(-i * H_Z)
    qml.ApproxTimeEvolution(H_Z, time=1, n=TROTTER_STEPS)   
    
    # Return state
    return qml.density_matrix(wires=[0, 1])

qfi_matrix = qml.gradients.quantum_fisher(circuit)(theta_values)
print(qfi_matrix)

Recap: I had to set diff_method=“parameter-shift” or 'finite-diff' to solve the

“shape-mismatch for sum”

issue related to using ApproxTimeEvolution and PhaseDamping simultaneously. However, when I set the hyperparameter n>1, I encounter the

“Cannot calculate ancestors for an operator that occurs multiple times.”