[BUG] `qml.eigvals` returns `nan` for a legal fractional power operator

[YujinSong-hep]

Expected behavior

qml.pow documents float exponents and explicitly shows fractional powers such as qml.pow(qml.X(0), 0.5).

For a legal fractional power operator, qml.eigvals(op) should return the eigenvalues of that operator. The result should agree with the eigenvalues of qml.matrix(op) up to ordering.

For example, qml.pow(qml.Z(0), 0.5) is the square root of a Pauli-Z phase operation. Since the eigenvalues of Z are 1 and -1, a valid square-root branch has finite complex eigenvalues 1 and 1j.

Actual behavior

qml.eigvals(qml.pow(qml.Z(0), 0.5)) returns a spectrum containing nan:

qml.eigvals: [np.float64(1.0), np.float64(nan)]
matrix eigvals: [1.000000e+00+0.j 6.123234e-17+1.j]

The same operator's matrix is well-defined and has finite complex eigenvalues.

Additional information

This appears to come from pennylane/ops/op_math/pow.py::Pow.eigvals, which computes

return [value**self.z for value in base_eigvals]

When value is a negative real eigenvalue and self.z is fractional, NumPy/Python real-power semantics produce nan instead of the expected complex root. The matrix path for the same operator already produces a finite complex operator matrix, so the eigenspectrum path is inconsistent with the operator itself.

Source code

import numpy as np
import pennylane as qml

op = qml.pow(qml.Z(0), 0.5)

print("qml.eigvals:", qml.eigvals(op))
print("matrix eigvals:", np.linalg.eigvals(qml.matrix(op)))

Tracebacks

Output：

qml.eigvals: [np.float64(1.0), np.float64(nan)]
matrix eigvals: [1.000000e+00+0.j 6.123234e-17+1.j]

System information

Version: 0.45.0
Platform info: macOS
Python version: 3.13.13

Existing GitHub issues

• I have searched existing GitHub issues to make sure the issue does not already exist.

[CatalinaAlbornoz]

Thanks for opening this bug report @YujinSong-hep !