diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 96613c882f1..dfb90a6acfb 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -2,6 +2,9 @@
 
 <h3>New features since last release</h3>
 
+* Added a new `SISWAP` operation and a `SQISW` alias with support to the `default_qubit` device.
+  [#1563](https://github.com/PennyLaneAI/pennylane/pull/1563)
+
 * The `RotosolveOptimizer` now can tackle general parametrized circuits, and is no longer
   restricted to single-qubit Pauli rotations.
   [(#1489)](https://github.com/PennyLaneAI/pennylane/pull/1489)
@@ -262,7 +265,7 @@ and requirements-ci.txt (unpinned). This latter would be used by the CI.
 This release contains contributions from (in alphabetical order):
 
 
-Akash Narayanan B, Thomas Bromley, Tanya Garg, Josh Izaac, Prateek Jain, Johannes Jakob Meyer, Maria Schuld,
+Vishnu Ajith, Akash Narayanan B, Thomas Bromley, Tanya Garg, Josh Izaac, Prateek Jain, Johannes Jakob Meyer, Maria Schuld,
 Ingrid Strandberg, David Wierichs, Vincent Wong.
 
 
diff --git a/doc/introduction/operations.rst b/doc/introduction/operations.rst
index 8a1a63894e8..0e08d343356 100644
--- a/doc/introduction/operations.rst
+++ b/doc/introduction/operations.rst
@@ -73,6 +73,8 @@ Qubit gates
     ~pennylane.CY
     ~pennylane.SWAP
     ~pennylane.ISWAP
+    ~pennylane.SISWAP
+    ~pennylane.SQISW
     ~pennylane.IsingXX
     ~pennylane.IsingYY
     ~pennylane.IsingZZ
diff --git a/pennylane/devices/default_qubit.py b/pennylane/devices/default_qubit.py
index ea8a60acfce..404e6c3d685 100644
--- a/pennylane/devices/default_qubit.py
+++ b/pennylane/devices/default_qubit.py
@@ -108,6 +108,8 @@ class DefaultQubit(QubitDevice):
         "CNOT",
         "SWAP",
         "ISWAP",
+        "SISWAP",
+        "SQISW",
         "CSWAP",
         "Toffoli",
         "CY",
diff --git a/pennylane/ops/qubit/__init__.py b/pennylane/ops/qubit/__init__.py
index e9a23e95ea2..729338fb184 100644
--- a/pennylane/ops/qubit/__init__.py
+++ b/pennylane/ops/qubit/__init__.py
@@ -52,6 +52,8 @@
     "CY",
     "SWAP",
     "ISWAP",
+    "SISWAP",
+    "SQISW",
     "CSWAP",
     "Toffoli",
     "RX",
diff --git a/pennylane/ops/qubit/non_parametric_ops.py b/pennylane/ops/qubit/non_parametric_ops.py
index 24fd38452ce..eb019773ed4 100644
--- a/pennylane/ops/qubit/non_parametric_ops.py
+++ b/pennylane/ops/qubit/non_parametric_ops.py
@@ -646,6 +646,69 @@ def adjoint(self):
         return ISWAP(wires=self.wires).inv()
 
 
+class SISWAP(Operation):
+    r"""SISWAP(wires)
+    The square root of i-swap operator. Can also be accessed as ``qml.SQISW``
+
+    .. math:: SISWAP = \begin{bmatrix}
+            1 & 0 & 0 & 0 \\
+            0 & 1/ \sqrt{2} & i/\sqrt{2} & 0\\
+            0 & i/ \sqrt{2} & 1/ \sqrt{2} & 0\\
+            0 & 0 & 0 & 1
+        \end{bmatrix}.
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 0
+
+    Args:
+        wires (Sequence[int]): the wires the operation acts on
+    """
+    num_params = 0
+    num_wires = 2
+    par_domain = None
+
+    @classmethod
+    def _matrix(cls, *params):
+        return np.array(
+            [
+                [1, 0, 0, 0],
+                [0, INV_SQRT2, INV_SQRT2 * 1j, 0],
+                [0, INV_SQRT2 * 1j, INV_SQRT2, 0],
+                [0, 0, 0, 1],
+            ]
+        )
+
+    @classmethod
+    def _eigvals(cls, *params):
+        return np.array([INV_SQRT2 * (1 + 1j), INV_SQRT2 * (1 - 1j), 1, 1])
+
+    @staticmethod
+    def decomposition(wires):
+        decomp_ops = [
+            SX(wires=wires[0]),
+            qml.RZ(np.pi / 2, wires=wires[0]),
+            CNOT(wires=[wires[0], wires[1]]),
+            SX(wires=wires[0]),
+            qml.RZ(7 * np.pi / 4, wires=wires[0]),
+            SX(wires=wires[0]),
+            qml.RZ(np.pi / 2, wires=wires[0]),
+            SX(wires=wires[1]),
+            qml.RZ(7 * np.pi / 4, wires=wires[1]),
+            CNOT(wires=[wires[0], wires[1]]),
+            SX(wires=wires[0]),
+            SX(wires=wires[1]),
+        ]
+        return decomp_ops
+
+    def adjoint(self):
+        return SISWAP(wires=self.wires).inv()
+
+
+SQISW = SISWAP
+
+
 class CSWAP(Operation):
     r"""CSWAP(wires)
     The controlled-swap operator
diff --git a/tests/devices/test_default_qubit.py b/tests/devices/test_default_qubit.py
index c5602efb12f..8498dc82d04 100644
--- a/tests/devices/test_default_qubit.py
+++ b/tests/devices/test_default_qubit.py
@@ -224,7 +224,57 @@ def test_apply_operation_single_wire_no_parameters_inverse(
         ),
     ]
 
-    all_two_wires_no_parameters = test_data_two_wires_no_parameters + test_data_iswap
+    test_data_siswap = [
+        (qml.SISWAP, [1, 0, 0, 0], [1, 0, 0, 0]),
+        (qml.SISWAP, [0, 1, 0, 0], [0, 1 / math.sqrt(2), 1 / math.sqrt(2) * 1j, 0]),
+        (
+            qml.SISWAP,
+            [1 / math.sqrt(2), 1 / math.sqrt(2), 0, 0],
+            [1 / math.sqrt(2), 0.5, 0.5 * 1j, 0],
+        ),
+    ]
+
+    test_data_sqisw = [
+        (qml.SQISW, [1, 0, 0, 0], [1, 0, 0, 0]),
+        (qml.SQISW, [0, 1, 0, 0], [0, 1 / math.sqrt(2), 1 / math.sqrt(2) * 1j, 0]),
+        (
+            qml.SQISW,
+            [1 / math.sqrt(2), 1 / math.sqrt(2), 0, 0],
+            [1 / math.sqrt(2), 0.5, 0.5 * 1j, 0],
+        ),
+    ]
+
+    test_data_siswap_inv = [
+        (
+            qml.SISWAP,
+            [1 / math.sqrt(2), 0, 1 / math.sqrt(2), 0],
+            [1 / math.sqrt(2), -0.5 * 1j, 0.5, 0],
+        ),
+        (qml.SISWAP, [0, 0, 1, 0], [0, -1 / math.sqrt(2) * 1j, 1 / math.sqrt(2), 0]),
+        (
+            qml.SISWAP,
+            [1 / math.sqrt(2), 0, -1 / math.sqrt(2), 0],
+            [1 / math.sqrt(2), 0.5 * 1j, -0.5, 0],
+        ),
+    ]
+
+    test_data_sqisw_inv = [
+        (
+            qml.SQISW,
+            [1 / math.sqrt(2), 0, 1 / math.sqrt(2), 0],
+            [1 / math.sqrt(2), -0.5 * 1j, 0.5, 0],
+        ),
+        (qml.SQISW, [0, 0, 1, 0], [0, -1 / math.sqrt(2) * 1j, 1 / math.sqrt(2), 0]),
+        (
+            qml.SQISW,
+            [1 / math.sqrt(2), 0, -1 / math.sqrt(2), 0],
+            [1 / math.sqrt(2), 0.5 * 1j, -0.5, 0],
+        ),
+    ]
+
+    all_two_wires_no_parameters = (
+        test_data_two_wires_no_parameters + test_data_iswap + test_data_siswap + test_data_sqisw
+    )
 
     @pytest.mark.parametrize("operation,input,expected_output", all_two_wires_no_parameters)
     def test_apply_operation_two_wires_no_parameters(
@@ -240,7 +290,12 @@ def test_apply_operation_two_wires_no_parameters(
             qubit_device_2_wires._state.flatten(), np.array(expected_output), atol=tol, rtol=0
         )
 
-    all_two_wires_no_parameters_inv = test_data_two_wires_no_parameters + test_data_iswap_inv
+    all_two_wires_no_parameters_inv = (
+        test_data_two_wires_no_parameters
+        + test_data_iswap_inv
+        + test_data_siswap_inv
+        + test_data_sqisw_inv
+    )
 
     @pytest.mark.parametrize("operation,input,expected_output", all_two_wires_no_parameters_inv)
     def test_apply_operation_two_wires_no_parameters_inverse(
diff --git a/tests/gate_data.py b/tests/gate_data.py
index dfae7c13ff4..9000f975742 100644
--- a/tests/gate_data.py
+++ b/tests/gate_data.py
@@ -23,6 +23,14 @@
 CNOT = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])  #: CNOT gate
 SWAP = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])  #: SWAP gate
 ISWAP = np.array([[1, 0, 0, 0], [0, 0, 1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]])  #: ISWAP gate
+SISWAP = np.array(
+    [
+        [1, 0, 0, 0],
+        [0, 1 / math.sqrt(2), 1 / math.sqrt(2) * 1j, 0],
+        [0, 1 / math.sqrt(2) * 1j, 1 / math.sqrt(2), 0],
+        [0, 0, 0, 1],
+    ]
+)
 CZ = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]])  #: CZ gate
 S = np.array([[1, 0], [0, 1j]])  #: Phase Gate
 T = np.array([[1, 0], [0, cmath.exp(1j * np.pi / 4)]])  #: T Gate
diff --git a/tests/ops/test_qubit_ops.py b/tests/ops/test_qubit_ops.py
index d71d4740057..5f2df5830f0 100644
--- a/tests/ops/test_qubit_ops.py
+++ b/tests/ops/test_qubit_ops.py
@@ -39,6 +39,7 @@
     CNOT,
     SWAP,
     ISWAP,
+    SISWAP,
     CZ,
     S,
     T,
@@ -400,6 +401,7 @@ def circuit(basis_state):
     (qml.CNOT, CNOT),
     (qml.SWAP, SWAP),
     (qml.ISWAP, ISWAP),
+    (qml.SISWAP, SISWAP),
     (qml.CZ, CZ),
     (qml.S, S),
     (qml.T, T),
@@ -419,6 +421,8 @@ def circuit(basis_state):
     qml.CY(wires=[0, 1]),
     qml.SWAP(wires=[0, 1]),
     qml.ISWAP(wires=[0, 1]),
+    qml.SISWAP(wires=[0, 1]),
+    qml.SQISW(wires=[0, 1]),
     qml.CSWAP(wires=[0, 1, 2]),
     qml.PauliRot(0.123, "Y", wires=0),
     qml.IsingXX(0.123, wires=[0, 1]),
@@ -772,6 +776,55 @@ def test_ISWAP_decomposition(self, tol):
 
         assert np.allclose(decomposed_matrix, op.matrix, atol=tol, rtol=0)
 
+    @pytest.mark.parametrize("siswap_op", [qml.SISWAP, qml.SQISW])
+    def test_SISWAP_decomposition(self, siswap_op, tol):
+        """Tests that the decomposition of the SISWAP gate and its SQISW alias gate is correct"""
+        op = siswap_op(wires=[0, 1])
+        res = op.decomposition(op.wires)
+
+        assert len(res) == 12
+
+        assert res[0].wires == Wires([0])
+        assert res[1].wires == Wires([0])
+        assert res[2].wires == Wires([0, 1])
+        assert res[3].wires == Wires([0])
+        assert res[4].wires == Wires([0])
+        assert res[5].wires == Wires([0])
+        assert res[6].wires == Wires([0])
+        assert res[7].wires == Wires([1])
+        assert res[8].wires == Wires([1])
+        assert res[9].wires == Wires([0, 1])
+        assert res[10].wires == Wires([0])
+        assert res[11].wires == Wires([1])
+
+        assert res[0].name == "SX"
+        assert res[1].name == "RZ"
+        assert res[2].name == "CNOT"
+        assert res[3].name == "SX"
+        assert res[4].name == "RZ"
+        assert res[5].name == "SX"
+        assert res[6].name == "RZ"
+        assert res[7].name == "SX"
+        assert res[8].name == "RZ"
+        assert res[9].name == "CNOT"
+        assert res[10].name == "SX"
+        assert res[11].name == "SX"
+
+        mats = []
+        for i in reversed(res):
+            if i.wires == Wires([1]):
+                mats.append(np.kron(np.eye(2), i.matrix))
+            elif i.wires == Wires([0]):
+                mats.append(np.kron(i.matrix, np.eye(2)))
+            elif i.wires == Wires([1, 0]) and i.name == "CNOT":
+                mats.append(np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]]))
+            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(self, tol):
         """Tests that the decomposition of the IsingXX gate is correct"""
         param = 0.1234
@@ -1536,6 +1589,14 @@ def test_iswap_eigenval(self):
         res = op.eigvals
         assert np.allclose(res, exp)
 
+    @pytest.mark.parametrize("siswap_op", [qml.SISWAP, qml.SQISW])
+    def test_siswap_eigenval(self, siswap_op):
+        """Tests that the ISWAP eigenvalue matches the numpy eigenvalues of the ISWAP matrix"""
+        op = siswap_op(wires=[0, 1])
+        exp = np.linalg.eigvals(op.matrix)
+        res = op.eigvals
+        assert np.allclose(res, exp)
+
     def test_swap_decomposition(self):
         """Tests the swap operator produces the correct output"""
         opr = qml.SWAP(wires=[0, 1])