diff --git a/fast_pauli/cpp/include/__summed_pauli_op.hpp b/fast_pauli/cpp/include/__summed_pauli_op.hpp
index 26bf0e4..856a9c0 100644
--- a/fast_pauli/cpp/include/__summed_pauli_op.hpp
+++ b/fast_pauli/cpp/include/__summed_pauli_op.hpp
@@ -531,6 +531,28 @@ template <std::floating_point T> struct SummedPauliOp
             }
         }
     }
+
+    /**
+     * @brief Get the dense representation of the SummedPauliOp as a 3D tensor
+     *
+     * @param A_k_out The output tensor to fill with the dense representation
+     */
+    void to_tensor(Tensor<3> A_k_out) const
+    {
+        for (size_t i = 0; i < pauli_strings.size(); ++i)
+        {
+            PauliString const &ps = pauli_strings[i];
+            auto [cols, vals] = get_sparse_repr<T>(ps.paulis);
+
+            for (size_t k = 0; k < n_operators(); ++k)
+            {
+                for (size_t j = 0; j < dim(); ++j)
+                {
+                    A_k_out(k, j, cols[j]) += coeffs(i, k) * vals[j];
+                }
+            }
+        }
+    }
 };
 
 } // namespace fast_pauli
diff --git a/fast_pauli/cpp/src/fast_pauli.cpp b/fast_pauli/cpp/src/fast_pauli.cpp
index 8cc3b85..a53cd36 100644
--- a/fast_pauli/cpp/src/fast_pauli.cpp
+++ b/fast_pauli/cpp/src/fast_pauli.cpp
@@ -1056,6 +1056,21 @@ Returns
 -------
 np.ndarray
     Expectation value(s) in a form of 2D numpy array (n_operators, n_states) according to the shape of input states
+)%")
+        .def(
+            "to_tensor",
+            [](fp::SummedPauliOp<float_type> const &self) {
+                auto dense_op =
+                    fp::__detail::owning_ndarray_from_shape<cfloat_t, 3>({self.n_operators(), self.dim(), self.dim()});
+                self.to_tensor(fp::__detail::ndarray_to_mdspan<cfloat_t, 3>(dense_op));
+                return dense_op;
+            },
+            R"%(Returns a dense representation of SummedPauliOp.
+
+Returns
+-------
+np.ndarray
+    3D numpy array of complex numbers with a shape of (n_operators, 2^n_qubits, 2^n_qubits)
 )%");
     //
     ;
diff --git a/tests/fast_pauli/test_summed_pauli_op.py b/tests/fast_pauli/test_summed_pauli_op.py
index d340164..7c81a38 100644
--- a/tests/fast_pauli/test_summed_pauli_op.py
+++ b/tests/fast_pauli/test_summed_pauli_op.py
@@ -95,3 +95,41 @@ def test_apply(
 
     # Check
     np.testing.assert_allclose(new_states, new_states_naive, atol=1e-13)
+
+
+@pytest.mark.parametrize(
+    "summed_pauli_op", [fp.SummedPauliOp], ids=resolve_parameter_repr
+)
+@pytest.mark.parametrize(
+    "n_states,n_operators,n_qubits",
+    [(s, o, q) for s in [1, 10, 1000] for o in [1, 10, 100] for q in [1, 2, 6]],
+)
+def test_to_tensor(
+    summed_pauli_op: type[fp.SummedPauliOp],
+    n_states: int,
+    n_operators: int,
+    n_qubits: int,
+) -> None:
+    """Test the dense representation of the SummedPauliOp."""
+    # initialize SummedPauliOp
+    pauli_strings = fp.helpers.calculate_pauli_strings_max_weight(n_qubits, 2)
+    n_strings = len(pauli_strings)
+    coeffs_2d = np.random.rand(n_strings, n_operators).astype(np.complex128)
+    op = summed_pauli_op(pauli_strings, coeffs_2d)
+
+    # get dense representation
+    dense_op = op.to_tensor()
+
+    # Check against doing it manually
+    # Get dense representation of each PauliString
+    ps_dense = np.array([ps.to_tensor() for ps in pauli_strings])
+
+    summed_pauli_op_check = np.zeros(
+        (n_operators, 2**n_qubits, 2**n_qubits), dtype=np.complex128
+    )
+
+    for k in range(n_operators):
+        for j in range(n_strings):
+            summed_pauli_op_check[k] += coeffs_2d[j, k] * ps_dense[j]
+
+    np.testing.assert_allclose(dense_op, summed_pauli_op_check)