diff --git a/doc/introduction/operations.rst b/doc/introduction/operations.rst
index 3facaa3f4ad..e1ac06ae033 100644
--- a/doc/introduction/operations.rst
+++ b/doc/introduction/operations.rst
@@ -155,6 +155,7 @@ Parametrized gates
     ~pennylane.U2
     ~pennylane.U3
     ~pennylane.IsingXX
+    ~pennylane.IsingXY
     ~pennylane.IsingYY
     ~pennylane.IsingZZ
 
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 8dbb86c9ce8..d6f6f5ea30a 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -298,6 +298,9 @@
 * Added separate requirements_dev.txt for separation of concerns for code development and just using PennyLane.
   [(#2635)](https://github.com/PennyLaneAI/pennylane/pull/2635)
 
+* Add `IsingXY` gate.
+  [(#2649)](https://github.com/PennyLaneAI/pennylane/pull/2649)
+
 * The performance of building sparse Hamiltonians has been improved by accumulating the sparse representation of coefficient-operator pairs in a temporary storage and by eliminating unnecessary `kron` operations on identity matrices. 
   [(#2630)](https://github.com/PennyLaneAI/pennylane/pull/2630)
 
@@ -392,6 +395,6 @@
 
 This release contains contributions from (in alphabetical order):
 
-Amintor Dusko, Chae-Yeun Park, Christian Gogolin, Christina Lee, David Wierichs, Edward Jiang, Guillermo Alonso-Linaje,
+Amintor Dusko, Ankit Khandelwal, Chae-Yeun Park, Christian Gogolin, Christina Lee, David Wierichs, Edward Jiang, Guillermo Alonso-Linaje,
 Jay Soni, Juan Miguel Arrazola, Katharine Hyatt, Korbinian, Kottmann, Maria Schuld, Mikhail Andrenkov, Romain Moyard,
 Qi Hu, Samuel Banning, Soran Jahangiri, Utkarsh Azad, WingCode
diff --git a/pennylane/devices/default_qubit.py b/pennylane/devices/default_qubit.py
index f3b3a1626dc..4a9b5cf0cd4 100644
--- a/pennylane/devices/default_qubit.py
+++ b/pennylane/devices/default_qubit.py
@@ -132,6 +132,7 @@ class DefaultQubit(QubitDevice):
         "IsingXX",
         "IsingYY",
         "IsingZZ",
+        "IsingXY",
         "SingleExcitation",
         "SingleExcitationPlus",
         "SingleExcitationMinus",
diff --git a/pennylane/devices/tests/test_gates.py b/pennylane/devices/tests/test_gates.py
index 2ede7210aa5..e4e42b87272 100755
--- a/pennylane/devices/tests/test_gates.py
+++ b/pennylane/devices/tests/test_gates.py
@@ -81,6 +81,7 @@
     "IsingXX": qml.IsingXX(0, wires=[0, 1]),
     "IsingYY": qml.IsingYY(0, wires=[0, 1]),
     "IsingZZ": qml.IsingZZ(0, wires=[0, 1]),
+    "IsingXY": qml.IsingXY(0, wires=[0, 1]),
     "SingleExcitation": qml.SingleExcitation(0, wires=[0, 1]),
     "SingleExcitationPlus": qml.SingleExcitationPlus(0, wires=[0, 1]),
     "SingleExcitationMinus": qml.SingleExcitationMinus(0, wires=[0, 1]),
@@ -192,6 +193,15 @@
     ]
 )
 
+IsingXY = lambda phi: np.array(
+    [
+        [1, 0, 0, 0],
+        [0, cos(phi / 2), 1j * sin(phi / 2), 0],
+        [0, 1j * sin(phi / 2), cos(phi / 2), 0],
+        [0, 0, 0, 1],
+    ]
+)
+
 IsingYY = lambda phi: np.array(
     [
         [cos(phi / 2), 0, 0, 1j * sin(phi / 2)],
@@ -255,6 +265,7 @@ def adjoint_tuple(op, orig_mat):
     (qml.CRY, cry),
     (qml.CRZ, crz),
     (qml.IsingXX, IsingXX),
+    (qml.IsingXY, IsingXY),
     (qml.IsingYY, IsingYY),
     (qml.IsingZZ, IsingZZ),
 ]
diff --git a/pennylane/ops/qubit/__init__.py b/pennylane/ops/qubit/__init__.py
index 404c0c351c8..9b0ae33508d 100644
--- a/pennylane/ops/qubit/__init__.py
+++ b/pennylane/ops/qubit/__init__.py
@@ -77,6 +77,7 @@
     "IsingXX",
     "IsingYY",
     "IsingZZ",
+    "IsingXY",
     "BasisState",
     "QubitStateVector",
     "QubitDensityMatrix",
diff --git a/pennylane/ops/qubit/parametric_ops.py b/pennylane/ops/qubit/parametric_ops.py
index 1c28c1f1443..20a59a0ab0e 100644
--- a/pennylane/ops/qubit/parametric_ops.py
+++ b/pennylane/ops/qubit/parametric_ops.py
@@ -2699,3 +2699,162 @@ def adjoint(self):
 
     def pow(self, z):
         return [IsingZZ(self.data[0] * z, wires=self.wires)]
+
+
+class IsingXY(Operation):
+    r"""
+    Ising XY coupling gate
+
+    .. math:: \mathtt{XY}(\phi) = \begin{bmatrix}
+            1 & 0 & 0 & 0 \\
+            0 & \cos(\phi / 2) & i \sin(\phi / 2) & 0 \\
+            0 & i \sin(\phi / 2) & \cos(\phi / 2) & 0 \\
+            0 & 0 & 0 & 1
+        \end{bmatrix}.
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 1
+    * Gradient recipe: The XY operator satisfies a four-term parameter-shift rule
+
+      .. math::
+          \frac{d}{d \phi} f(XY(\phi))
+          = c_+ \left[ f(XY(\phi + a)) - f(XY(\phi - a)) \right]
+          - c_- \left[ f(XY(\phi + b)) - f(XY(\phi - b)) \right]
+
+      where :math:`f` is an expectation value depending on :math:`XY(\phi)`, and
+
+      - :math:`a = \pi / 2`
+      - :math:`b = 3 \pi / 2`
+      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4 \sqrt{2}}`
+
+    Args:
+        phi (float): the phase angle
+        wires (int): the subsystem the gate acts on
+        do_queue (bool): Indicates whether the operator should be
+            immediately pushed into the Operator queue (optional)
+        id (str or None): String representing the operation (optional)
+    """
+    num_wires = 2
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    grad_method = "A"
+    parameter_frequencies = [(0.5, 1.0)]
+
+    def generator(self):
+        return 0.25 * qml.PauliX(wires=self.wires[0]) @ qml.PauliX(
+            wires=self.wires[1]
+        ) + 0.25 * qml.PauliY(wires=self.wires[0]) @ qml.PauliY(wires=self.wires[1])
+
+    def __init__(self, phi, wires, do_queue=True, id=None):
+        super().__init__(phi, wires=wires, do_queue=do_queue, id=id)
+
+    @staticmethod
+    def compute_decomposition(phi, wires):
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+
+        .. seealso:: :meth:`~.IsingXY.decomposition`.
+
+        Args:
+            phi (float): the phase angle
+            wires (Iterable, Wires): the subsystem the gate acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.IsingXY.compute_decomposition(1.23, wires=(0,1))
+        [Hadamard(wires=[0]), CY(wires=[0, 1]), RY(0.615, wires=[0]), RX(-0.615, wires=[1]), CY(wires=[0, 1]), Hadamard(wires=[0])]
+
+        """
+        return [
+            qml.Hadamard(wires=[wires[0]]),
+            qml.CY(wires=wires),
+            qml.RY(phi / 2, wires=[wires[0]]),
+            qml.RX(-phi / 2, wires=[wires[1]]),
+            qml.CY(wires=wires),
+            qml.Hadamard(wires=[wires[0]]),
+        ]
+
+    @staticmethod
+    def compute_matrix(phi):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.IsingXY.matrix`
+
+
+        Args:
+           phi (tensor_like or float): phase angle
+
+        Returns:
+           tensor_like: canonical matrix
+
+        **Example**
+
+        >>> qml.IsingXY.compute_matrix(0.5)
+        array([[1.        +0.j        , 0.        +0.j        ,        0.        +0.j        , 0.        +0.j        ],
+               [0.        +0.j        , 0.96891242+0.j        ,        0.        +0.24740396j, 0.        +0.j        ],
+               [0.        +0.j        , 0.        +0.24740396j,        0.96891242+0.j        , 0.        +0.j        ],
+               [0.        +0.j        , 0.        +0.j        ,        0.        +0.j        , 1.        +0.j        ]])
+        """
+        c = qml.math.cos(phi / 2)
+        s = qml.math.sin(phi / 2)
+        Y = qml.math.convert_like(np.diag([0, 1, 1, 0])[::-1].copy(), phi)
+
+        if qml.math.get_interface(phi) == "tensorflow":
+            c = qml.math.cast_like(c, 1j)
+            s = qml.math.cast_like(s, 1j)
+            Y = qml.math.cast_like(Y, 1j)
+
+        return qml.math.diag([1, c, c, 1]) + 1j * s * Y
+
+    @staticmethod
+    def compute_eigvals(phi):  # pylint: disable=arguments-differ
+        r"""Eigenvalues of the operator in the computational basis (static method).
+
+        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U`,
+        the operator can be reconstructed as
+
+        .. math:: O = U \Sigma U^{\dagger},
+
+        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
+
+        Otherwise, no particular order for the eigenvalues is guaranteed.
+
+        .. seealso:: :meth:`~.IsingXY.eigvals`
+
+
+        Args:
+            phi (tensor_like or float): phase angle
+
+        Returns:
+            tensor_like: eigenvalues
+
+        **Example**
+
+        >>> qml.IsingXY.compute_eigvals(0.5)
+        array([0.96891242+0.24740396j, 0.96891242-0.24740396j,       1.        +0.j        , 1.        +0.j        ])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            phi = qml.math.cast_like(phi, 1j)
+
+        pos_phase = qml.math.exp(1.0j * phi / 2)
+        neg_phase = qml.math.exp(-1.0j * phi / 2)
+
+        return qml.math.stack([pos_phase, neg_phase, 1, 1])
+
+    def adjoint(self):
+        (phi,) = self.parameters
+        return IsingXY(-phi, wires=self.wires)
+
+    def pow(self, z):
+        return [IsingXY(self.data[0] * z, wires=self.wires)]
diff --git a/tests/gate_data.py b/tests/gate_data.py
index 76707f593ed..190327c0a22 100644
--- a/tests/gate_data.py
+++ b/tests/gate_data.py
@@ -279,6 +279,29 @@ def IsingYY(phi):
     return np.cos(phi / 2) * II - 1j * np.sin(phi / 2) * YY
 
 
+def IsingXY(phi):
+    r"""Ising XY coupling gate.
+
+    .. math:: \mathtt{XY}(\phi) = \begin{bmatrix}
+            1 & 0 & 0 & 0 \\
+            0 & \cos(\phi / 2) & i \sin(\phi / 2) & 0 \\
+            0 & i \sin(\phi / 2) & \cos(\phi / 2) & 0 \\
+            0 & 0 & 0 & 1
+        \end{bmatrix}.
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+    Returns:
+        array[complex]: unitary 4x4 rotation matrix
+    """
+    mat = II.copy()
+    mat[1][1] = np.cos(phi / 2)
+    mat[2][2] = np.cos(phi / 2)
+    mat[1][2] = 1j * np.sin(phi / 2)
+    mat[2][1] = 1j * np.sin(phi / 2)
+    return mat
+
+
 def IsingZZ(phi):
     r"""Ising ZZ coupling gate
 
diff --git a/tests/ops/qubit/test_parametric_ops.py b/tests/ops/qubit/test_parametric_ops.py
index 2ab2670889e..1a950d68f9b 100644
--- a/tests/ops/qubit/test_parametric_ops.py
+++ b/tests/ops/qubit/test_parametric_ops.py
@@ -33,6 +33,7 @@
     qml.IsingXX(0.123, wires=[0, 1]),
     qml.IsingYY(0.123, wires=[0, 1]),
     qml.IsingZZ(0.123, wires=[0, 1]),
+    qml.IsingXY(0.123, wires=[0, 1]),
     qml.Rot(0.123, 0.456, 0.789, wires=0),
     qml.PhaseShift(2.133, wires=0),
     qml.ControlledPhaseShift(1.777, wires=[0, 2]),
@@ -526,6 +527,43 @@ def test_isingxx_decomposition(self, tol):
 
         assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
 
+    def test_isingxy_decomposition(self, tol):
+        """Tests that the decomposition of the IsingXY gate is correct"""
+        param = 0.1234
+        op = qml.IsingXY(param, wires=[3, 2])
+        res = op.decomposition()
+
+        assert len(res) == 6
+
+        assert res[0].wires == Wires([3])
+        assert res[1].wires == Wires([3, 2])
+        assert res[2].wires == Wires([3])
+        assert res[3].wires == Wires([2])
+        assert res[4].wires == Wires([3, 2])
+        assert res[5].wires == Wires([3])
+
+        assert res[0].name == "Hadamard"
+        assert res[1].name == "CY"
+        assert res[2].name == "RY"
+        assert res[3].name == "RX"
+        assert res[4].name == "CY"
+        assert res[5].name == "Hadamard"
+
+        mats = []
+        for i in reversed(res):
+            if i.wires == Wires([3]):
+                # RY and Hadamard gate
+                mats.append(np.kron(i.matrix(), np.eye(2)))
+            elif i.wires == Wires([2]):
+                # RX gate
+                mats.append(np.kron(np.eye(2), i.matrix()))
+            else:
+                mats.append(i.matrix())
+
+        decomposed_matrix = np.linalg.multi_dot(mats)
+
+        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
+
     def test_isingxx_decomposition_broadcasted(self, tol):
         """Tests that the decomposition of the broadcasted IsingXX gate is correct"""
         param = np.array([-0.1, 0.2, 0.5])
@@ -881,6 +919,91 @@ def get_expected(theta):
             qml.IsingXX(param, wires=[0, 1]).matrix(), get_expected(param), atol=tol, rtol=0
         )
 
+    def test_isingxy(self, tol):
+        """Test that the IsingXY operation is correct"""
+        assert np.allclose(qml.IsingXY.compute_matrix(0), np.identity(4), atol=tol, rtol=0)
+        assert np.allclose(qml.IsingXY(0, wires=[0, 1]).matrix(), np.identity(4), atol=tol, rtol=0)
+
+        def get_expected(theta):
+            expected = np.eye(4, dtype=np.complex128)
+            expected[1][1] = np.cos(theta / 2)
+            expected[2][2] = np.cos(theta / 2)
+            expected[1][2] = 1j * np.sin(theta / 2)
+            expected[2][1] = 1j * np.sin(theta / 2)
+            return expected
+
+        param = np.pi / 2
+        assert np.allclose(qml.IsingXY.compute_matrix(param), get_expected(param), atol=tol, rtol=0)
+        assert np.allclose(
+            qml.IsingXY(param, wires=[0, 1]).matrix(), get_expected(param), atol=tol, rtol=0
+        )
+
+        param = np.pi
+        assert np.allclose(qml.IsingXY.compute_matrix(param), get_expected(param), atol=tol, rtol=0)
+        assert np.allclose(
+            qml.IsingXY(param, wires=[0, 1]).matrix(), get_expected(param), atol=tol, rtol=0
+        )
+
+    @pytest.mark.parametrize("phi", np.linspace(-np.pi, np.pi, 10))
+    def test_isingxy_eigvals(self, phi, tol):
+        """Test eigenvalues computation for IsingXY"""
+        evs = qml.IsingXY.compute_eigvals(phi)
+        evs_expected = [
+            qml.math.cos(phi / 2) + 1j * qml.math.sin(phi / 2),
+            qml.math.cos(phi / 2) - 1j * qml.math.sin(phi / 2),
+            1,
+            1,
+        ]
+        assert qml.math.allclose(evs, evs_expected)
+
+    @pytest.mark.tf
+    @pytest.mark.parametrize("phi", np.linspace(-np.pi, np.pi, 10))
+    def test_isingxy_eigvals_tf(self, phi, tol):
+        """Test eigenvalues computation for IsingXY using Tensorflow interface"""
+        import tensorflow as tf
+
+        param_tf = tf.Variable(phi)
+        evs = qml.IsingXY.compute_eigvals(param_tf)
+        evs_expected = [
+            qml.math.cos(phi / 2) + 1j * qml.math.sin(phi / 2),
+            qml.math.cos(phi / 2) - 1j * qml.math.sin(phi / 2),
+            1,
+            1,
+        ]
+        assert qml.math.allclose(evs, evs_expected)
+
+    @pytest.mark.torch
+    @pytest.mark.parametrize("phi", np.linspace(-np.pi, np.pi, 10))
+    def test_isingxy_eigvals_torch(self, phi, tol):
+        """Test eigenvalues computation for IsingXY using Torch interface"""
+        import torch
+
+        param_torch = torch.tensor(phi)
+        evs = qml.IsingXY.compute_eigvals(param_torch)
+        evs_expected = [
+            qml.math.cos(phi / 2) + 1j * qml.math.sin(phi / 2),
+            qml.math.cos(phi / 2) - 1j * qml.math.sin(phi / 2),
+            1,
+            1,
+        ]
+        assert qml.math.allclose(evs, evs_expected)
+
+    @pytest.mark.jax
+    @pytest.mark.parametrize("phi", np.linspace(-np.pi, np.pi, 10))
+    def test_isingxy_eigvals_jax(self, phi, tol):
+        """Test eigenvalues computation for IsingXY using JAX interface"""
+        import jax
+
+        param_jax = jax.numpy.array(phi)
+        evs = qml.IsingXY.compute_eigvals(param_jax)
+        evs_expected = [
+            qml.math.cos(phi / 2) + 1j * qml.math.sin(phi / 2),
+            qml.math.cos(phi / 2) - 1j * qml.math.sin(phi / 2),
+            1,
+            1,
+        ]
+        assert qml.math.allclose(evs, evs_expected)
+
     def test_isingxx_broadcasted(self, tol):
         """Test that the broadcasted IsingXX operation is correct"""
         z = np.zeros(3)
@@ -1520,6 +1643,72 @@ def circuit(phi):
         res = qml.grad(circuit)(phi)
         assert np.allclose(res, expected, atol=tol, rtol=0)
 
+    @pytest.mark.autograd
+    @pytest.mark.parametrize("dev_name,diff_method,phi", configuration)
+    def test_isingxy_autograd_grad(self, tol, dev_name, diff_method, phi):
+        """Test the gradient with Autograd for the gate IsingXY."""
+        dev = qml.device(dev_name, wires=2)
+
+        psi_0 = 0.1
+        psi_1 = 0.2
+        psi_2 = 0.3
+        psi_3 = 0.4
+
+        init_state = npp.array([psi_0, psi_1, psi_2, psi_3], requires_grad=False)
+        norm = np.linalg.norm(init_state)
+        init_state /= norm
+
+        @qml.qnode(dev, diff_method=diff_method, interface="autograd")
+        def circuit(phi):
+            qml.QubitStateVector(init_state, wires=[0, 1])
+            qml.IsingXY(phi, wires=[0, 1])
+            return qml.expval(qml.PauliZ(0))
+
+        phi = npp.array(0.1, requires_grad=True)
+
+        expected = (1 / norm**2) * (psi_2**2 - psi_1**2) * np.sin(phi)
+
+        res = qml.grad(circuit)(phi)
+        assert np.allclose(res, expected, atol=tol, rtol=0)
+
+    @pytest.mark.jax
+    @pytest.mark.parametrize("dev_name,diff_method,phi", configuration)
+    def test_isingxy_jax_grad(self, tol, dev_name, diff_method, phi):
+        """Test the gradient with JAX for the gate IsingXY."""
+
+        if diff_method in {"finite-diff"}:
+            pytest.skip("Test does not support finite-diff")
+
+        if diff_method in {"parameter-shift"}:
+            pytest.skip("Test does not support parameter-shift")
+
+        import jax
+        import jax.numpy as jnp
+
+        dev = qml.device(dev_name, wires=2)
+
+        psi_0 = 0.1
+        psi_1 = 0.2
+        psi_2 = 0.3
+        psi_3 = 0.4
+
+        init_state = jnp.array([psi_0, psi_1, psi_2, psi_3])
+        norm = jnp.linalg.norm(init_state)
+        init_state = init_state / norm
+
+        @qml.qnode(dev, diff_method=diff_method, interface="jax")
+        def circuit(phi):
+            qml.QubitStateVector(init_state, wires=[0, 1])
+            qml.IsingXY(phi, wires=[0, 1])
+            return qml.expval(qml.PauliZ(0))
+
+        phi = jnp.array(0.1)
+
+        expected = (1 / norm**2) * (psi_2**2 - psi_1**2) * np.sin(phi)
+
+        res = jax.grad(circuit, argnums=0)(phi)
+        assert np.allclose(res, expected, atol=tol, rtol=0)
+
     @pytest.mark.jax
     @pytest.mark.parametrize("dev_name,diff_method,phi", configuration)
     def test_isingxx_jax_grad(self, tol, dev_name, diff_method, phi):
@@ -1654,6 +1843,38 @@ def circuit(phi):
         res = jax.grad(circuit, argnums=0)(phi)
         assert np.allclose(res, expected, atol=tol, rtol=0)
 
+    @pytest.mark.tf
+    @pytest.mark.parametrize("dev_name,diff_method,phi", configuration)
+    def test_isingxy_tf_grad(self, tol, dev_name, diff_method, phi):
+        """Test the gradient with Tensorflow for the gate IsingXY."""
+        import tensorflow as tf
+
+        dev = qml.device(dev_name, wires=2)
+
+        psi_0 = tf.Variable(0.1, dtype=tf.complex128)
+        psi_1 = tf.Variable(0.2, dtype=tf.complex128)
+        psi_2 = tf.Variable(0.3, dtype=tf.complex128)
+        psi_3 = tf.Variable(0.4, dtype=tf.complex128)
+
+        init_state = tf.Variable([psi_0, psi_1, psi_2, psi_3], dtype=tf.complex128)
+        norm = tf.norm(init_state)
+        init_state = init_state / norm
+
+        @qml.qnode(dev, interface="tf", diff_method=diff_method)
+        def circuit(phi):
+            qml.QubitStateVector(init_state, wires=[0, 1])
+            qml.IsingXY(phi, wires=[0, 1])
+            return qml.expval(qml.PauliZ(0))
+
+        phi = tf.Variable(0.1, dtype=tf.complex128)
+
+        expected = (1 / norm**2) * (psi_2**2 - psi_1**2) * tf.sin(phi)
+
+        with tf.GradientTape() as tape:
+            result = circuit(phi)
+        res = tape.gradient(result, phi)
+        assert np.allclose(res, expected, atol=tol, rtol=0)
+
     @pytest.mark.tf
     @pytest.mark.parametrize("dev_name,diff_method,phi", configuration)
     def test_isingxx_tf_grad(self, tol, dev_name, diff_method, phi):
@@ -2606,6 +2827,7 @@ def test_string_parameter_broadcasted(self):
     qml.U1(1.23, wires=0),
     qml.IsingXX(-2.345, wires=(0, 1)),
     qml.IsingYY(3.1652, wires=(0, 1)),
+    qml.IsingXY(-1.234, wires=(0, 1)),
     qml.IsingZZ(1.789, wires=("a", "b")),
     # broadcasted ops
     qml.RX(np.array([1.234, 4.129]), wires=0),
@@ -2659,6 +2881,8 @@ def test_pow_matrix(self, op, n):
     (qml.U2(1.234, 2.345, wires=0), Wires([])),
     (qml.U3(1.234, 2.345, 3.456, wires=0), Wires([])),
     (qml.IsingXX(1.234, wires=(0, 1)), Wires([])),
+    (qml.IsingYY(1.234, wires=(0, 1)), Wires([])),
+    (qml.IsingXY(1.234, wires=(0, 1)), Wires([])),
     (qml.IsingYY(np.array([-5.1, 0.219]), wires=(0, 1)), Wires([])),
     (qml.IsingZZ(1.234, wires=(0, 1)), Wires([])),
     ### Controlled Ops