diff --git a/CMakeLists.txt b/CMakeLists.txt
index b61d7fba9..01867014a 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -63,6 +63,29 @@ option(MQT_CORE_INSTALL "Generate installation instructions for MQT Core"
 option(BUILD_MQT_CORE_TESTS "Build tests for the MQT Core project" ${MQT_CORE_MASTER_PROJECT})
 option(BUILD_MQT_CORE_SHARED_LIBS "Build MQT Core libraries as shared libraries"
        ${BUILD_SHARED_LIBS})
+option(BUILD_MQT_CORE_MLIR "Build the MLIR submodule of the MQT Core project" ON)
+if(APPLE AND CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
+  set(BUILD_MQT_CORE_MLIR
+      OFF
+      CACHE BOOL
+            "Disable MLIR build on macOS with GCC to avoid ABI issues between STL and libstdc++"
+            FORCE)
+  message(
+    WARNING "Disabling MLIR build on macOS with GCC to avoid ABI issues between STL and libstdc++")
+endif()
+if(APPLE AND CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
+  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "17.0.0")
+    set(BUILD_MQT_CORE_MLIR
+        OFF
+        CACHE
+          BOOL
+          "Disable MLIR build on macOS with Apple Clang < 17 due to missing complete C++20 support"
+          FORCE)
+    message(
+      WARNING
+        "Disabling MLIR build on macOS with Apple Clang < 17 due to missing complete C++20 support")
+  endif()
+endif()
 
 # try to determine the project version
 include(GetVersion)
@@ -93,29 +116,6 @@ endif()
 
 set(MQT_CORE_TARGET_NAME "mqt-core")
 
-option(BUILD_MQT_CORE_MLIR "Build the MLIR submodule of the MQT Core project" ON)
-if(APPLE AND CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
-  set(BUILD_MQT_CORE_MLIR
-      OFF
-      CACHE BOOL
-            "Disable MLIR build on macOS with GCC to avoid ABI issues between STL and libstdc++"
-            FORCE)
-  message(
-    WARNING "Disabling MLIR build on macOS with GCC to avoid ABI issues between STL and libstdc++")
-endif()
-if(APPLE AND CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
-  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "17.0.0")
-    set(BUILD_MQT_CORE_MLIR
-        OFF
-        CACHE
-          BOOL
-          "Disable MLIR build on macOS with Apple Clang < 17 due to missing complete C++20 support"
-          FORCE)
-    message(
-      WARNING
-        "Disabling MLIR build on macOS with Apple Clang < 17 due to missing complete C++20 support")
-  endif()
-endif()
 if(BUILD_MQT_CORE_MLIR)
   include(SetupMLIR)
 endif()
diff --git a/cmake/ExternalDependencies.cmake b/cmake/ExternalDependencies.cmake
index 1176e8fa4..7f352efb3 100644
--- a/cmake/ExternalDependencies.cmake
+++ b/cmake/ExternalDependencies.cmake
@@ -22,6 +22,19 @@ if(BUILD_MQT_CORE_BINDINGS)
   find_package(nanobind CONFIG REQUIRED)
 endif()
 
+if(BUILD_MQT_CORE_MLIR)
+  set(Eigen_VERSION
+      5.0.1
+      CACHE STRING "Eigen version")
+  set(Eigen_URL
+      https://gitlab.com/libeigen/eigen/-/archive/${Eigen_VERSION}/eigen-${Eigen_VERSION}.tar.gz)
+  set(EIGEN_BUILD_TESTING
+      OFF
+      CACHE INTERNAL "Disable building Eigen tests")
+  FetchContent_Declare(Eigen URL ${Eigen_URL} FIND_PACKAGE_ARGS ${Eigen_VERSION})
+  list(APPEND FETCH_PACKAGES Eigen)
+endif()
+
 set(JSON_VERSION
     3.12.0
     CACHE STRING "nlohmann_json version")
diff --git a/mlir/.clang-tidy b/mlir/.clang-tidy
index 9f2230b31..68b1fc7bc 100644
--- a/mlir/.clang-tidy
+++ b/mlir/.clang-tidy
@@ -7,3 +7,7 @@ Checks: |
   -misc-use-anonymous-namespace,
   llvm-twine-local,
   -cppcoreguidelines-pro-bounds-avoid-unchecked-container-access
+
+CheckOptions:
+  - key: misc-include-cleaner.IgnoreHeaders
+    value: "Eigen/.*;unsupported/Eigen/.*"
diff --git a/mlir/include/mlir/Dialect/QCO/IR/QCODialect.h b/mlir/include/mlir/Dialect/QCO/IR/QCODialect.h
index 11fdbb5a6..67fd0c12a 100644
--- a/mlir/include/mlir/Dialect/QCO/IR/QCODialect.h
+++ b/mlir/include/mlir/Dialect/QCO/IR/QCODialect.h
@@ -21,6 +21,7 @@
 #pragma clang diagnostic pop
 #endif
 
+#include <Eigen/Core>
 #include <mlir/Bytecode/BytecodeOpInterface.h>
 #include <mlir/Dialect/Arith/IR/Arith.h>
 #include <mlir/IR/Value.h>
diff --git a/mlir/include/mlir/Dialect/QCO/IR/QCOInterfaces.td b/mlir/include/mlir/Dialect/QCO/IR/QCOInterfaces.td
index 5ae789c1d..b7e08461f 100644
--- a/mlir/include/mlir/Dialect/QCO/IR/QCOInterfaces.td
+++ b/mlir/include/mlir/Dialect/QCO/IR/QCOInterfaces.td
@@ -28,6 +28,62 @@ def UnitaryOpInterface : OpInterface<"UnitaryOpInterface"> {
 
     let cppNamespace = "::mlir::qco";
 
+    // Generic implementation body for getUnitaryMatrix methods
+    defvar unitaryMatrixMethodBody = [{
+      auto process = [&]<typename MatrixType>(MatrixType&& m) -> bool {
+          using TargetT = std::remove_cvref_t<decltype(out)>;
+          using SourceT = std::remove_cvref_t<MatrixType>;
+
+          constexpr bool isTargetDynamic =
+              (TargetT::SizeAtCompileTime == Eigen::Dynamic);
+          constexpr bool isSourceDynamic =
+              (SourceT::SizeAtCompileTime == Eigen::Dynamic);
+
+          // Case 1: Target is Dynamic. Always accepts source.
+          if constexpr (isTargetDynamic) {
+            out = std::forward<MatrixType>(m);
+            return true;
+          }
+          // Case 2: Target is Fixed.
+          else {
+            // Case 2a: Source is Dynamic. Runtime dimension check required.
+            if constexpr (isSourceDynamic) {
+              if (m.rows() == static_cast<Eigen::Index>(TargetT::RowsAtCompileTime) &&
+                  m.cols() == static_cast<Eigen::Index>(TargetT::ColsAtCompileTime))
+                  [[likely]] {
+                out = std::forward<MatrixType>(m);
+                return true;
+              }
+            }
+            // Case 2b: Source is Fixed. Compile-time check.
+            else if constexpr (static_cast<Eigen::Index>(
+                                   SourceT::RowsAtCompileTime) ==
+                                   static_cast<Eigen::Index>(
+                                       TargetT::RowsAtCompileTime) &&
+                               static_cast<Eigen::Index>(
+                                   SourceT::ColsAtCompileTime) ==
+                                   static_cast<Eigen::Index>(
+                                       TargetT::ColsAtCompileTime)) {
+              out = std::forward<MatrixType>(m);
+              return true;
+            }
+          }
+          return false;
+        };
+
+
+      if constexpr (requires { $_op.getUnitaryMatrix().has_value(); }) {
+        if (auto&& matrix = $_op.getUnitaryMatrix()) {
+          return process(std::move(*matrix));
+        }
+        return false;
+      } else if constexpr (requires { $_op.getUnitaryMatrix(); }) {
+        return process($_op.getUnitaryMatrix());
+      } else {
+        llvm::reportFatalUsageError("Operation '" + $_op.getBaseSymbol() + "' has no unitary matrix definition!");
+      }
+    }];
+
     let methods = [
         // Qubit accessors
         InterfaceMethod<
@@ -107,7 +163,69 @@ def UnitaryOpInterface : OpInterface<"UnitaryOpInterface"> {
             "Returns the base symbol/mnemonic of the operation.",
             "StringRef", "getBaseSymbol", (ins)
         >,
+        InterfaceMethod<
+            "Populates the given 1x1 unitary matrix if possible.",
+            "bool", "getUnitaryMatrix1x1",
+            (ins "Eigen::Matrix<std::complex<double>, 1, 1>&":$out),
+            unitaryMatrixMethodBody
+        >,
+        InterfaceMethod<
+            "Populates the given 2x2 unitary matrix if possible.",
+            "bool", "getUnitaryMatrix2x2",
+            (ins "Eigen::Matrix2cd&":$out),
+            unitaryMatrixMethodBody
+        >,
+        InterfaceMethod<
+            "Populates the given 4x4 unitary matrix if possible.",
+            "bool", "getUnitaryMatrix4x4",
+            (ins "Eigen::Matrix4cd&":$out),
+            unitaryMatrixMethodBody
+        >,
+        InterfaceMethod<
+            "Populates the given dynamic unitary matrix.",
+            "bool", "getUnitaryMatrixDynamic",
+            (ins "Eigen::MatrixXcd&":$out),
+            unitaryMatrixMethodBody
+        >
     ];
+
+    let extraClassDeclaration = [{
+        template<typename MatrixType>
+        std::optional<MatrixType> getUnitaryMatrix() {
+          MatrixType out;
+          bool result = false;
+
+          // Dispatch to the appropriate fixed-size or dynamic method based on the
+          // matrix type.
+          if constexpr (MatrixType::RowsAtCompileTime == 1 &&
+                        MatrixType::ColsAtCompileTime == 1) {
+            result = this->getUnitaryMatrix1x1(out);
+          } else if constexpr (MatrixType::RowsAtCompileTime == 2 &&
+                               MatrixType::ColsAtCompileTime == 2) {
+            result = this->getUnitaryMatrix2x2(out);
+          } else if constexpr (MatrixType::RowsAtCompileTime == 4 &&
+                               MatrixType::ColsAtCompileTime == 4) {
+            result = this->getUnitaryMatrix4x4(out);
+          } else if constexpr (MatrixType::SizeAtCompileTime == Eigen::Dynamic) {
+            result = this->getUnitaryMatrixDynamic(out);
+          } else {
+            // Fallback: Try obtaining dynamic matrix and see if size matches
+            Eigen::MatrixXcd dynamicOut;
+            if (this->getUnitaryMatrixDynamic(dynamicOut)) {
+              if (dynamicOut.rows() == MatrixType::RowsAtCompileTime &&
+                  dynamicOut.cols() == MatrixType::ColsAtCompileTime) {
+                out = dynamicOut;
+                result = true;
+              }
+            }
+          }
+
+          if (result) {
+            return out;
+          }
+          return std::nullopt;
+        }
+    }];
 }
 
 #endif // QCO_INTERFACES
diff --git a/mlir/include/mlir/Dialect/QCO/IR/QCOOps.td b/mlir/include/mlir/Dialect/QCO/IR/QCOOps.td
index aef50a319..1ccf16320 100644
--- a/mlir/include/mlir/Dialect/QCO/IR/QCOOps.td
+++ b/mlir/include/mlir/Dialect/QCO/IR/QCOOps.td
@@ -239,7 +239,7 @@ def ResetOp : QCOOp<"reset", [Idempotent, SameOperandsAndResultType]> {
 //===----------------------------------------------------------------------===//
 
 class TargetAndParameterArityTrait<int T, int P>
-        : ParamNativeOpTrait<"TargetAndParameterArityTrait", !strconcat(!cast<string>(T), ",", !cast<string>(P))> {
+        : ParamNativeOpTrait<"TargetAndParameterArityTrait", !cast<string>(T) # ", " # !cast<string>(P)> {
     let cppNamespace = "::mlir::qco";
 }
 
@@ -271,7 +271,8 @@ def GPhaseOp : QCOOp<"gphase", traits = [UnitaryOpInterface, ZeroTargetOneParame
     let assemblyFormat = "`(` $theta `)` attr-dict";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "gphase"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "gphase"; }
+        [[nodiscard]] std::optional<Eigen::Matrix<std::complex<double>, 1, 1>> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -297,7 +298,8 @@ def IdOp : QCOOp<"id", traits = [UnitaryOpInterface, OneTargetZeroParameter]> {
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "id"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "id"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -319,7 +321,8 @@ def XOp : QCOOp<"x", traits = [UnitaryOpInterface, OneTargetZeroParameter]> {
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "x"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "x"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -341,7 +344,8 @@ def YOp : QCOOp<"y", traits = [UnitaryOpInterface, OneTargetZeroParameter]> {
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "y"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "y"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -363,7 +367,8 @@ def ZOp : QCOOp<"z", traits = [UnitaryOpInterface, OneTargetZeroParameter]> {
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "z"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "z"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -385,7 +390,8 @@ def HOp : QCOOp<"h", traits = [UnitaryOpInterface, OneTargetZeroParameter]> {
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "h"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "h"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -407,7 +413,8 @@ def SOp : QCOOp<"s", traits = [UnitaryOpInterface, OneTargetZeroParameter]> {
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "s"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "s"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -429,7 +436,8 @@ def SdgOp : QCOOp<"sdg", traits = [UnitaryOpInterface, OneTargetZeroParameter]>
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "sdg"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "sdg"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -451,7 +459,8 @@ def TOp : QCOOp<"t", traits = [UnitaryOpInterface, OneTargetZeroParameter]> {
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "t"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "t"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -473,7 +482,8 @@ def TdgOp : QCOOp<"tdg", traits = [UnitaryOpInterface, OneTargetZeroParameter]>
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "tdg"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "tdg"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -495,7 +505,8 @@ def SXOp : QCOOp<"sx", traits = [UnitaryOpInterface, OneTargetZeroParameter]> {
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "sx"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "sx"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -517,7 +528,8 @@ def SXdgOp : QCOOp<"sxdg", traits = [UnitaryOpInterface, OneTargetZeroParameter]
     let assemblyFormat = "$qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "sxdg"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "sxdg"; }
+        [[nodiscard]] static Eigen::Matrix2cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -540,7 +552,8 @@ def RXOp : QCOOp<"rx", traits = [UnitaryOpInterface, OneTargetOneParameter]> {
     let assemblyFormat = "`(` $theta `)` $qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "rx"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "rx"; }
+        [[nodiscard]] std::optional<Eigen::Matrix2cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -567,7 +580,8 @@ def RYOp : QCOOp<"ry", traits = [UnitaryOpInterface, OneTargetOneParameter]> {
     let assemblyFormat = "`(` $theta `)` $qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "ry"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "ry"; }
+        [[nodiscard]] std::optional<Eigen::Matrix2cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -594,7 +608,8 @@ def RZOp : QCOOp<"rz", traits = [UnitaryOpInterface, OneTargetOneParameter]> {
     let assemblyFormat = "`(` $theta `)` $qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "rz"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "rz"; }
+        [[nodiscard]] std::optional<Eigen::Matrix2cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -621,7 +636,8 @@ def POp : QCOOp<"p", traits = [UnitaryOpInterface, OneTargetOneParameter]> {
     let assemblyFormat = "`(` $theta `)` $qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "p"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "p"; }
+        [[nodiscard]] std::optional<Eigen::Matrix2cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -649,7 +665,8 @@ def ROp : QCOOp<"r", traits = [UnitaryOpInterface, OneTargetTwoParameter]> {
     let assemblyFormat = "`(` $theta `,` $phi `)` $qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "r"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "r"; }
+        [[nodiscard]] std::optional<Eigen::Matrix2cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -677,7 +694,8 @@ def U2Op : QCOOp<"u2", traits = [UnitaryOpInterface, OneTargetTwoParameter]> {
     let assemblyFormat = "`(` $phi `,` $lambda `)` $qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "u2"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "u2"; }
+        [[nodiscard]] std::optional<Eigen::Matrix2cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -706,7 +724,8 @@ def UOp : QCOOp<"u", traits = [UnitaryOpInterface, OneTargetThreeParameter]> {
     let assemblyFormat = "`(` $theta `,` $phi `,` $lambda `)` $qubit_in attr-dict `:` type($qubit_in) `->` type($qubit_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "u"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "u"; }
+        [[nodiscard]] std::optional<Eigen::Matrix2cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -733,7 +752,8 @@ def SWAPOp : QCOOp<"swap", traits = [UnitaryOpInterface, TwoTargetZeroParameter]
     let assemblyFormat = "$qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "swap"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "swap"; }
+        [[nodiscard]] static Eigen::Matrix4cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -756,7 +776,8 @@ def iSWAPOp : QCOOp<"iswap", traits = [UnitaryOpInterface, TwoTargetZeroParamete
     let assemblyFormat = "$qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "iswap"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "iswap"; }
+        [[nodiscard]] static Eigen::Matrix4cd getUnitaryMatrix();
     }];
 }
 
@@ -777,7 +798,8 @@ def DCXOp : QCOOp<"dcx", traits = [UnitaryOpInterface, TwoTargetZeroParameter]>
     let assemblyFormat = "$qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "dcx"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "dcx"; }
+        [[nodiscard]] static Eigen::Matrix4cd getUnitaryMatrix();
     }];
 }
 
@@ -798,7 +820,8 @@ def ECROp : QCOOp<"ecr", traits = [UnitaryOpInterface, TwoTargetZeroParameter]>
     let assemblyFormat = "$qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "ecr"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "ecr"; }
+        [[nodiscard]] static Eigen::Matrix4cd getUnitaryMatrix();
     }];
 
     let hasCanonicalizer = 1;
@@ -822,7 +845,8 @@ def RXXOp : QCOOp<"rxx", traits = [UnitaryOpInterface, TwoTargetOneParameter]> {
     let assemblyFormat = "`(` $theta `)` $qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "rxx"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "rxx"; }
+        [[nodiscard]] std::optional<Eigen::Matrix4cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -850,7 +874,8 @@ def RYYOp : QCOOp<"ryy", traits = [UnitaryOpInterface, TwoTargetOneParameter]> {
     let assemblyFormat = "`(` $theta `)` $qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "ryy"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "ryy"; }
+        [[nodiscard]] std::optional<Eigen::Matrix4cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -878,7 +903,8 @@ def RZXOp : QCOOp<"rzx", traits = [UnitaryOpInterface, TwoTargetOneParameter]> {
     let assemblyFormat = "`(` $theta `)` $qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "rzx"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "rzx"; }
+        [[nodiscard]] std::optional<Eigen::Matrix4cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -906,7 +932,8 @@ def RZZOp : QCOOp<"rzz", traits = [UnitaryOpInterface, TwoTargetOneParameter]> {
     let assemblyFormat = "`(` $theta `)` $qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "rzz"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "rzz"; }
+        [[nodiscard]] std::optional<Eigen::Matrix4cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -935,7 +962,8 @@ def XXPlusYYOp : QCOOp<"xx_plus_yy", traits = [UnitaryOpInterface, TwoTargetTwoP
     let assemblyFormat = "`(` $theta `,` $beta `)` $qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "xx_plus_yy"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "xx_plus_yy"; }
+        [[nodiscard]] std::optional<Eigen::Matrix4cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -964,7 +992,8 @@ def XXMinusYYOp : QCOOp<"xx_minus_yy", traits = [UnitaryOpInterface, TwoTargetTw
     let assemblyFormat = "`(` $theta `,` $beta `)` $qubit0_in `,` $qubit1_in attr-dict `:` type($qubit0_in) `,` type($qubit1_in) `->` type($qubit0_out) `,` type($qubit1_out)";
 
     let extraClassDeclaration = [{
-        static StringRef getBaseSymbol() { return "xx_minus_yy"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "xx_minus_yy"; }
+        [[nodiscard]] std::optional<Eigen::Matrix4cd> getUnitaryMatrix();
     }];
 
     let builders = [
@@ -1003,7 +1032,7 @@ def BarrierOp : QCOOp<"barrier", traits = [UnitaryOpInterface]> {
         Value getOutputForInput(Value input);
         static size_t getNumParams();
         static Value getParameter(size_t i);
-        static StringRef getBaseSymbol() { return "barrier"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "barrier"; }
     }];
 
     let builders = [
@@ -1093,7 +1122,8 @@ def CtrlOp : QCOOp<"ctrl", traits =
         Value getOutputForInput(Value input);
         size_t getNumParams();
         Value getParameter(size_t i);
-        static StringRef getBaseSymbol() { return "ctrl"; }
+        [[nodiscard]] static StringRef getBaseSymbol() { return "ctrl"; }
+        [[nodiscard]] std::optional<Eigen::MatrixXcd> getUnitaryMatrix();
     }];
 
     let builders = [
diff --git a/mlir/lib/Dialect/QCO/IR/CMakeLists.txt b/mlir/lib/Dialect/QCO/IR/CMakeLists.txt
index c41cffd76..7c62ed47d 100644
--- a/mlir/lib/Dialect/QCO/IR/CMakeLists.txt
+++ b/mlir/lib/Dialect/QCO/IR/CMakeLists.txt
@@ -27,6 +27,7 @@ add_mlir_dialect_library(
   MLIRArithDialect
   MLIRInferTypeOpInterface
   MLIRSideEffectInterfaces
+  Eigen3::Eigen
   DISABLE_INSTALL)
 
 # collect header files
diff --git a/mlir/lib/Dialect/QCO/IR/Modifiers/CtrlOp.cpp b/mlir/lib/Dialect/QCO/IR/Modifiers/CtrlOp.cpp
index 5d49b3a44..caa07a6e7 100644
--- a/mlir/lib/Dialect/QCO/IR/Modifiers/CtrlOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Modifiers/CtrlOp.cpp
@@ -10,7 +10,9 @@
 
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 
+#include <cassert>
 #include <cstddef>
+#include <cstdint>
 #include <llvm/ADT/STLExtras.h>
 #include <llvm/ADT/STLFunctionalExtras.h>
 #include <llvm/ADT/SmallVector.h>
@@ -24,6 +26,7 @@
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LLVM.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 
 using namespace mlir;
 using namespace mlir::qco;
@@ -363,3 +366,30 @@ void CtrlOp::getCanonicalizationPatterns(RewritePatternSet& results,
       .add<MergeNestedCtrl, RemoveTrivialCtrl, CtrlInlineGPhase, CtrlInlineId>(
           context);
 }
+
+std::optional<Eigen::MatrixXcd> CtrlOp::getUnitaryMatrix() {
+  auto&& bodyUnitary = getBodyUnitary();
+  if (!bodyUnitary) {
+    return std::nullopt;
+  }
+  auto&& targetMatrix = bodyUnitary.getUnitaryMatrix<Eigen::MatrixXcd>();
+  if (!targetMatrix) {
+    return std::nullopt;
+  }
+
+  // get dimensions of target matrix
+  const auto targetDim = targetMatrix->cols();
+  assert(targetMatrix->cols() == targetMatrix->rows());
+
+  // define dimensions and type of output matrix
+  assert(getNumControls() < sizeof(unsigned long long) * 8);
+  const auto dim = static_cast<int64_t>((1ULL << getNumControls()) * targetDim);
+
+  // initialize result with identity
+  Eigen::MatrixXcd matrix = Eigen::MatrixXcd::Identity(dim, dim);
+
+  // apply target matrix
+  matrix.bottomRightCorner(targetDim, targetDim) = *targetMatrix;
+
+  return matrix;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/DCXOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/DCXOp.cpp
new file mode 100644
index 000000000..7267c1de0
--- /dev/null
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/DCXOp.cpp
@@ -0,0 +1,23 @@
+/*
+ * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
+ * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
+ * All rights reserved.
+ *
+ * SPDX-License-Identifier: MIT
+ *
+ * Licensed under the MIT License
+ */
+
+#include "mlir/Dialect/QCO/IR/QCODialect.h"
+
+#include <Eigen/Core>
+
+using namespace mlir;
+using namespace mlir::qco;
+
+Eigen::Matrix4cd DCXOp::getUnitaryMatrix() {
+  return Eigen::Matrix4cd{{1, 0, 0, 0},  // row 0
+                          {0, 0, 1, 0},  // row 1
+                          {0, 0, 0, 1},  // row 2
+                          {0, 1, 0, 0}}; // row 3
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ECROp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ECROp.cpp
index fb204bd24..d7c9728a2 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ECROp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ECROp.cpp
@@ -11,10 +11,13 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <numbers>
 
 using namespace mlir;
 using namespace mlir::qco;
@@ -39,3 +42,15 @@ void ECROp::getCanonicalizationPatterns(RewritePatternSet& results,
                                         MLIRContext* context) {
   results.add<RemoveSubsequentECR>(context);
 }
+
+Eigen::Matrix4cd ECROp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  constexpr auto m0 = 0i;
+  constexpr auto m1 = std::complex<double>{1.0 / std::numbers::sqrt2};
+  constexpr auto mi = std::complex<double>{0.0, 1.0 / std::numbers::sqrt2};
+  return Eigen::Matrix4cd{{m0, m0, m1, mi},   // row 0
+                          {m0, m0, mi, m1},   // row 1
+                          {m1, -mi, m0, m0},  // row 2
+                          {-mi, m1, m0, m0}}; // row 3
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/GPhaseOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/GPhaseOp.cpp
index 677658c84..b23f88ef5 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/GPhaseOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/GPhaseOp.cpp
@@ -11,12 +11,15 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
 #include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -33,8 +36,8 @@ struct RemoveTrivialGPhase final : OpRewritePattern<GPhaseOp> {
 
   LogicalResult matchAndRewrite(GPhaseOp op,
                                 PatternRewriter& rewriter) const override {
-    const auto theta = valueToDouble(op.getTheta());
-    if (!theta || std::abs(*theta) > TOLERANCE) {
+    if (const auto theta = valueToDouble(op.getTheta());
+        !theta || std::abs(*theta) > TOLERANCE) {
       return failure();
     }
 
@@ -47,7 +50,8 @@ struct RemoveTrivialGPhase final : OpRewritePattern<GPhaseOp> {
 
 void GPhaseOp::build(OpBuilder& odsBuilder, OperationState& odsState,
                      const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, thetaOperand);
 }
 
@@ -55,3 +59,11 @@ void GPhaseOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                            MLIRContext* context) {
   results.add<RemoveTrivialGPhase>(context);
 }
+
+std::optional<Eigen::Matrix<std::complex<double>, 1, 1>>
+GPhaseOp::getUnitaryMatrix() {
+  if (const auto theta = valueToDouble(getTheta())) {
+    return Eigen::Matrix<std::complex<double>, 1, 1>{std::polar(1.0, *theta)};
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/HOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/HOp.cpp
index 30eeabb67..ca1f9ab53 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/HOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/HOp.cpp
@@ -11,10 +11,12 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <numbers>
 
 using namespace mlir;
 using namespace mlir::qco;
@@ -39,3 +41,8 @@ void HOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<RemoveSubsequentH>(context);
 }
+
+Eigen::Matrix2cd HOp::getUnitaryMatrix() {
+  constexpr auto x = 1.0 / std::numbers::sqrt2;
+  return Eigen::Matrix2cd{{x, x}, {x, -1.0 * x}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/IdOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/IdOp.cpp
index 555b05213..c1c1274c8 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/IdOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/IdOp.cpp
@@ -10,6 +10,7 @@
 
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 
+#include <Eigen/Core>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -39,3 +40,7 @@ void IdOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                        MLIRContext* context) {
   results.add<RemoveId>(context);
 }
+
+Eigen::Matrix2cd IdOp::getUnitaryMatrix() {
+  return Eigen::Matrix2cd{{1, 0}, {0, 1}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/POp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/POp.cpp
index 21f1ac0eb..9628bd373 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/POp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/POp.cpp
@@ -12,11 +12,14 @@
 #include "mlir/Dialect/QCO/QCOUtils.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -54,7 +57,8 @@ struct RemoveTrivialP final : OpRewritePattern<POp> {
 
 void POp::build(OpBuilder& odsBuilder, OperationState& odsState, Value qubitIn,
                 const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, qubitIn, thetaOperand);
 }
 
@@ -62,3 +66,10 @@ void POp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<MergeSubsequentP, RemoveTrivialP>(context);
 }
+
+std::optional<Eigen::Matrix2cd> POp::getUnitaryMatrix() {
+  if (const auto theta = valueToDouble(getTheta())) {
+    return Eigen::Matrix2cd{{1.0, 0.0}, {0.0, std::polar(1.0, *theta)}};
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ROp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ROp.cpp
index 56c7fc4c1..3eab2953a 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ROp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ROp.cpp
@@ -11,13 +11,16 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
 #include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
 #include <numbers>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -73,8 +76,9 @@ struct ReplaceRWithRY final : OpRewritePattern<ROp> {
 void ROp::build(OpBuilder& odsBuilder, OperationState& odsState, Value qubitIn,
                 const std::variant<double, Value>& theta,
                 const std::variant<double, Value>& phi) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
-  auto phiOperand = variantToValue(odsBuilder, odsState.location, phi);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
+  const auto phiOperand = variantToValue(odsBuilder, odsState.location, phi);
   build(odsBuilder, odsState, qubitIn, thetaOperand, phiOperand);
 }
 
@@ -82,3 +86,17 @@ void ROp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<ReplaceRWithRX, ReplaceRWithRY>(context);
 }
+
+std::optional<Eigen::Matrix2cd> ROp::getUnitaryMatrix() {
+  const auto theta = valueToDouble(getTheta());
+  const auto phi = valueToDouble(getPhi());
+  if (!theta || !phi) {
+    return std::nullopt;
+  }
+
+  const auto thetaSin = std::sin(*theta / 2.0);
+  const auto m01 = std::polar(thetaSin, -*phi - (std::numbers::pi / 2));
+  const auto m10 = std::polar(thetaSin, *phi - (std::numbers::pi / 2));
+  const std::complex<double> thetaCos = std::cos(*theta / 2.0);
+  return Eigen::Matrix2cd{{thetaCos, m01}, {m10, thetaCos}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RXOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RXOp.cpp
index 8abff26b3..b8f0aa3c5 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RXOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RXOp.cpp
@@ -12,11 +12,15 @@
 #include "mlir/Dialect/QCO/QCOUtils.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
+#include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -54,7 +58,8 @@ struct RemoveTrivialRX final : OpRewritePattern<RXOp> {
 
 void RXOp::build(OpBuilder& odsBuilder, OperationState& odsState, Value qubitIn,
                  const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, qubitIn, thetaOperand);
 }
 
@@ -62,3 +67,14 @@ void RXOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                        MLIRContext* context) {
   results.add<MergeSubsequentRX, RemoveTrivialRX>(context);
 }
+
+std::optional<Eigen::Matrix2cd> RXOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  if (const auto theta = valueToDouble(getTheta())) {
+    const auto m00 = std::cos(*theta / 2.0) + 0i;
+    const auto m01 = -1i * std::sin(*theta / 2.0);
+    return Eigen::Matrix2cd{{m00, m01}, {m01, m00}};
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RXXOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RXXOp.cpp
index ed2b74109..88f257083 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RXXOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RXXOp.cpp
@@ -12,11 +12,15 @@
 #include "mlir/Dialect/QCO/QCOUtils.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
+#include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -55,7 +59,8 @@ struct RemoveTrivialRXX final : OpRewritePattern<RXXOp> {
 void RXXOp::build(OpBuilder& odsBuilder, OperationState& odsState,
                   Value qubit0In, Value qubit1In,
                   const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, qubit0In, qubit1In, thetaOperand);
 }
 
@@ -63,3 +68,18 @@ void RXXOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                         MLIRContext* context) {
   results.add<MergeSubsequentRXX, RemoveTrivialRXX>(context);
 }
+
+std::optional<Eigen::Matrix4cd> RXXOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  if (const auto theta = valueToDouble(getTheta())) {
+    const auto m0 = 0i;
+    const auto mc = std::cos(*theta / 2.0) + 0i;
+    const auto ms = -1i * std::sin(*theta / 2.0);
+    return Eigen::Matrix4cd{{mc, m0, m0, ms},  // row 0
+                            {m0, mc, ms, m0},  // row 1
+                            {m0, ms, mc, m0},  // row 2
+                            {ms, m0, m0, mc}}; // row 3
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RYOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RYOp.cpp
index d3037ca70..c4b7ef194 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RYOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RYOp.cpp
@@ -12,11 +12,15 @@
 #include "mlir/Dialect/QCO/QCOUtils.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
+#include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -54,7 +58,8 @@ struct RemoveTrivialRY final : OpRewritePattern<RYOp> {
 
 void RYOp::build(OpBuilder& odsBuilder, OperationState& odsState, Value qubitIn,
                  const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, qubitIn, thetaOperand);
 }
 
@@ -62,3 +67,12 @@ void RYOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                        MLIRContext* context) {
   results.add<MergeSubsequentRY, RemoveTrivialRY>(context);
 }
+
+std::optional<Eigen::Matrix2cd> RYOp::getUnitaryMatrix() {
+  if (const auto theta = valueToDouble(getTheta())) {
+    const auto m00 = std::complex<double>{std::cos(*theta / 2.0)};
+    const auto m01 = std::complex<double>{-std::sin(*theta / 2.0)};
+    return Eigen::Matrix2cd{{m00, m01}, {-m01, m00}};
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RYYOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RYYOp.cpp
index b30b5ae86..0a0940189 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RYYOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RYYOp.cpp
@@ -12,11 +12,15 @@
 #include "mlir/Dialect/QCO/QCOUtils.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
+#include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -55,7 +59,8 @@ struct RemoveTrivialRYY final : OpRewritePattern<RYYOp> {
 void RYYOp::build(OpBuilder& odsBuilder, OperationState& odsState,
                   Value qubit0In, Value qubit1In,
                   const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, qubit0In, qubit1In, thetaOperand);
 }
 
@@ -63,3 +68,18 @@ void RYYOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                         MLIRContext* context) {
   results.add<MergeSubsequentRYY, RemoveTrivialRYY>(context);
 }
+
+std::optional<Eigen::Matrix4cd> RYYOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  if (const auto theta = valueToDouble(getTheta())) {
+    const auto m0 = 0i;
+    const auto mc = std::complex<double>{std::cos(*theta / 2.0)};
+    const auto ms = std::complex<double>{0.0, std::sin(*theta / 2.0)};
+    return Eigen::Matrix4cd{{mc, m0, m0, ms},  // row 0
+                            {m0, mc, -ms, m0}, // row 1
+                            {m0, -ms, mc, m0}, // row 2
+                            {ms, m0, m0, mc}}; // row 3
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZOp.cpp
index 282f16d72..030647afc 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZOp.cpp
@@ -12,11 +12,14 @@
 #include "mlir/Dialect/QCO/QCOUtils.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -54,7 +57,8 @@ struct RemoveTrivialRZ final : OpRewritePattern<RZOp> {
 
 void RZOp::build(OpBuilder& odsBuilder, OperationState& odsState, Value qubitIn,
                  const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, qubitIn, thetaOperand);
 }
 
@@ -62,3 +66,15 @@ void RZOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                        MLIRContext* context) {
   results.add<MergeSubsequentRZ, RemoveTrivialRZ>(context);
 }
+
+std::optional<Eigen::Matrix2cd> RZOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  if (const auto theta = valueToDouble(getTheta())) {
+    const auto m00 = std::polar(1.0, -*theta / 2.0);
+    const auto m01 = 0i;
+    const auto m11 = std::polar(1.0, *theta / 2.0);
+    return Eigen::Matrix2cd{{m00, m01}, {m01, m11}};
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZXOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZXOp.cpp
index ab1f948ed..289b51110 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZXOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZXOp.cpp
@@ -12,11 +12,15 @@
 #include "mlir/Dialect/QCO/QCOUtils.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
+#include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -55,7 +59,8 @@ struct RemoveTrivialRZX final : OpRewritePattern<RZXOp> {
 void RZXOp::build(OpBuilder& odsBuilder, OperationState& odsState,
                   Value qubit0In, Value qubit1In,
                   const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, qubit0In, qubit1In, thetaOperand);
 }
 
@@ -63,3 +68,18 @@ void RZXOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                         MLIRContext* context) {
   results.add<MergeSubsequentRZX, RemoveTrivialRZX>(context);
 }
+
+std::optional<Eigen::Matrix4cd> RZXOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  if (const auto theta = valueToDouble(getTheta())) {
+    const auto m0 = 0i;
+    const auto mc = std::complex<double>{std::cos(*theta / 2.0)};
+    const auto ms = std::complex<double>{0.0, std::sin(*theta / 2.0)};
+    return Eigen::Matrix4cd{{mc, -ms, m0, m0}, // row 0
+                            {-ms, mc, m0, m0}, // row 1
+                            {m0, m0, mc, ms},  // row 2
+                            {m0, m0, ms, mc}}; // row 3
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZZOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZZOp.cpp
index 90f74f798..f16dbf74a 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZZOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/RZZOp.cpp
@@ -12,11 +12,14 @@
 #include "mlir/Dialect/QCO/QCOUtils.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -55,7 +58,8 @@ struct RemoveTrivialRZZ final : OpRewritePattern<RZZOp> {
 void RZZOp::build(OpBuilder& odsBuilder, OperationState& odsState,
                   Value qubit0In, Value qubit1In,
                   const std::variant<double, Value>& theta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
   build(odsBuilder, odsState, qubit0In, qubit1In, thetaOperand);
 }
 
@@ -63,3 +67,18 @@ void RZZOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                         MLIRContext* context) {
   results.add<MergeSubsequentRZZ, RemoveTrivialRZZ>(context);
 }
+
+std::optional<Eigen::Matrix4cd> RZZOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  if (const auto theta = valueToDouble(getTheta())) {
+    const auto m0 = 0i;
+    const auto mp = std::polar(1.0, *theta / 2.0);
+    const auto mm = std::polar(1.0, -*theta / 2.0);
+    return Eigen::Matrix4cd{{mm, m0, m0, m0},  // row 0
+                            {m0, mp, m0, m0},  // row 1
+                            {m0, m0, mp, m0},  // row 2
+                            {m0, m0, m0, mm}}; // row 3
+  }
+  return std::nullopt;
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SOp.cpp
index e3a435b53..623756878 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SOp.cpp
@@ -11,6 +11,8 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -51,3 +53,10 @@ void SOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<RemoveSAfterSdg, MergeSubsequentS>(context);
 }
+
+Eigen::Matrix2cd SOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  return Eigen::Matrix2cd{{1.0, 0.0}, // row 0
+                          {0.0, 1i}}; // row 1
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SWAPOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SWAPOp.cpp
index ffcd32c13..5a9b0bf24 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SWAPOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SWAPOp.cpp
@@ -11,6 +11,7 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -39,3 +40,10 @@ void SWAPOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                          MLIRContext* context) {
   results.add<RemoveSubsequentSWAP>(context);
 }
+
+Eigen::Matrix4cd SWAPOp::getUnitaryMatrix() {
+  return Eigen::Matrix4cd{{1, 0, 0, 0},  // row 0
+                          {0, 0, 1, 0},  // row 1
+                          {0, 1, 0, 0},  // row 2
+                          {0, 0, 0, 1}}; // row 3
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SXOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SXOp.cpp
index 465e01aff..8a0608ef7 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SXOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SXOp.cpp
@@ -11,6 +11,8 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -51,3 +53,9 @@ void SXOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                        MLIRContext* context) {
   results.add<RemoveSXAfterSXdg, MergeSubsequentSX>(context);
 }
+
+Eigen::Matrix2cd SXOp::getUnitaryMatrix() {
+  constexpr auto m00 = std::complex<double>{0.5, 0.5};
+  constexpr auto m01 = std::complex<double>{0.5, -0.5};
+  return Eigen::Matrix2cd{{m00, m01}, {m01, m00}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SXdgOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SXdgOp.cpp
index 73dc6ad8a..1136e090f 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SXdgOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SXdgOp.cpp
@@ -11,6 +11,8 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -52,3 +54,9 @@ void SXdgOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                          MLIRContext* context) {
   results.add<RemoveSXdgAfterSX, MergeSubsequentSXdg>(context);
 }
+
+Eigen::Matrix2cd SXdgOp::getUnitaryMatrix() {
+  constexpr auto m00 = std::complex<double>{0.5, -0.5};
+  constexpr auto m01 = std::complex<double>{0.5, 0.5};
+  return Eigen::Matrix2cd{{m00, m01}, {m01, m00}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SdgOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SdgOp.cpp
index 68ad4dce8..1c2c61425 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SdgOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/SdgOp.cpp
@@ -11,6 +11,8 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -51,3 +53,10 @@ void SdgOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                         MLIRContext* context) {
   results.add<RemoveSdgAfterS, MergeSubsequentSdg>(context);
 }
+
+Eigen::Matrix2cd SdgOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  return Eigen::Matrix2cd{{1.0, 0.0},  // row 0
+                          {0.0, -1i}}; // row 1
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/TOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/TOp.cpp
index c0bd0c7dc..32d70f445 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/TOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/TOp.cpp
@@ -11,10 +11,13 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <numbers>
 
 using namespace mlir;
 using namespace mlir::qco;
@@ -51,3 +54,8 @@ void TOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<RemoveTAfterTdg, MergeSubsequentT>(context);
 }
+
+Eigen::Matrix2cd TOp::getUnitaryMatrix() {
+  const auto m11 = std::polar(1.0, std::numbers::pi / 4.0);
+  return Eigen::Matrix2cd{{1.0, 0.0}, {0.0, m11}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/TdgOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/TdgOp.cpp
index dca978187..378dc0000 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/TdgOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/TdgOp.cpp
@@ -11,10 +11,13 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <numbers>
 
 using namespace mlir;
 using namespace mlir::qco;
@@ -52,3 +55,8 @@ void TdgOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                         MLIRContext* context) {
   results.add<RemoveTdgAfterT, MergeSubsequentTdg>(context);
 }
+
+Eigen::Matrix2cd TdgOp::getUnitaryMatrix() {
+  const auto m11 = std::polar(1.0, -std::numbers::pi / 4.0);
+  return Eigen::Matrix2cd{{1.0, 0.0}, {0.0, m11}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/U2Op.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/U2Op.cpp
index f9ef3ad37..6872bc559 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/U2Op.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/U2Op.cpp
@@ -11,13 +11,16 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
 #include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
 #include <numbers>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -99,8 +102,9 @@ struct ReplaceU2WithRY final : OpRewritePattern<U2Op> {
 void U2Op::build(OpBuilder& odsBuilder, OperationState& odsState, Value qubitIn,
                  const std::variant<double, Value>& phi,
                  const std::variant<double, Value>& lambda) {
-  auto phiOperand = variantToValue(odsBuilder, odsState.location, phi);
-  auto lambdaOperand = variantToValue(odsBuilder, odsState.location, lambda);
+  const auto phiOperand = variantToValue(odsBuilder, odsState.location, phi);
+  const auto lambdaOperand =
+      variantToValue(odsBuilder, odsState.location, lambda);
   build(odsBuilder, odsState, qubitIn, phiOperand, lambdaOperand);
 }
 
@@ -108,3 +112,20 @@ void U2Op::getCanonicalizationPatterns(RewritePatternSet& results,
                                        MLIRContext* context) {
   results.add<ReplaceU2WithH, ReplaceU2WithRX, ReplaceU2WithRY>(context);
 }
+
+std::optional<Eigen::Matrix2cd> U2Op::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  const auto phi = valueToDouble(getPhi());
+  const auto lambda = valueToDouble(getLambda());
+  if (!phi || !lambda) {
+    return std::nullopt;
+  }
+
+  const auto m00 = 1.0 / std::numbers::sqrt2 + 0i;
+  const auto m01 =
+      std::polar(1.0 / std::numbers::sqrt2, *lambda + std::numbers::pi);
+  const auto m10 = std::polar(1.0 / std::numbers::sqrt2, *phi);
+  const auto m11 = std::polar(1.0 / std::numbers::sqrt2, *phi + *lambda);
+  return Eigen::Matrix2cd{{m00, m01}, {m10, m11}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/UOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/UOp.cpp
index b4e783c20..37c0aa49a 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/UOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/UOp.cpp
@@ -11,13 +11,16 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
 #include <cmath>
+#include <complex>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
 #include <numbers>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -101,9 +104,11 @@ void UOp::build(OpBuilder& odsBuilder, OperationState& odsState, Value qubitIn,
                 const std::variant<double, Value>& theta,
                 const std::variant<double, Value>& phi,
                 const std::variant<double, Value>& lambda) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
-  auto phiOperand = variantToValue(odsBuilder, odsState.location, phi);
-  auto lambdaOperand = variantToValue(odsBuilder, odsState.location, lambda);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
+  const auto phiOperand = variantToValue(odsBuilder, odsState.location, phi);
+  const auto lambdaOperand =
+      variantToValue(odsBuilder, odsState.location, lambda);
   build(odsBuilder, odsState, qubitIn, thetaOperand, phiOperand, lambdaOperand);
 }
 
@@ -111,3 +116,22 @@ void UOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<ReplaceUWithP, ReplaceUWithRX, ReplaceUWithRY>(context);
 }
+
+std::optional<Eigen::Matrix2cd> UOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  const auto theta = valueToDouble(getTheta());
+  const auto phi = valueToDouble(getPhi());
+  const auto lambda = valueToDouble(getLambda());
+  if (!theta || !phi || !lambda) {
+    return std::nullopt;
+  }
+
+  const auto c = std::cos(*theta / 2.0);
+  const auto s = std::sin(*theta / 2.0);
+  const auto m00 = c + 0i;
+  const auto m01 = std::polar(s, *lambda + std::numbers::pi);
+  const auto m10 = std::polar(s, *phi);
+  const auto m11 = std::polar(c, *phi + *lambda);
+  return Eigen::Matrix2cd{{m00, m01}, {m10, m11}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XOp.cpp
index eb40915cb..e19075f78 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XOp.cpp
@@ -11,6 +11,7 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -39,3 +40,7 @@ void XOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<RemoveSubsequentX>(context);
 }
+
+Eigen::Matrix2cd XOp::getUnitaryMatrix() {
+  return Eigen::Matrix2cd{{0, 1}, {1, 0}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XXMinusYYOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XXMinusYYOp.cpp
index 9adba3254..be89e8895 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XXMinusYYOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XXMinusYYOp.cpp
@@ -11,7 +11,9 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
 #include <cmath>
+#include <complex>
 #include <mlir/Dialect/Arith/IR/Arith.h>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/BuiltinAttributes.h>
@@ -19,6 +21,8 @@
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <numbers>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -75,8 +79,9 @@ void XXMinusYYOp::build(OpBuilder& odsBuilder, OperationState& odsState,
                         Value qubit0In, Value qubit1In,
                         const std::variant<double, Value>& theta,
                         const std::variant<double, Value>& beta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
-  auto betaOperand = variantToValue(odsBuilder, odsState.location, beta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
+  const auto betaOperand = variantToValue(odsBuilder, odsState.location, beta);
   build(odsBuilder, odsState, qubit0In, qubit1In, thetaOperand, betaOperand);
 }
 
@@ -84,3 +89,24 @@ void XXMinusYYOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                               MLIRContext* context) {
   results.add<MergeSubsequentXXMinusYY>(context);
 }
+
+std::optional<Eigen::Matrix4cd> XXMinusYYOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  const auto theta = valueToDouble(getTheta());
+  const auto beta = valueToDouble(getBeta());
+  if (!theta || !beta) {
+    return std::nullopt;
+  }
+
+  const auto m0 = 0.0 + 0i;
+  const auto m1 = 1.0 + 0i;
+  const auto mc = std::cos(*theta / 2.0) + 0i;
+  const auto s = std::sin(*theta / 2.0);
+  const auto msp = std::polar(s, *beta - (std::numbers::pi / 2.));
+  const auto msm = std::polar(s, -*beta - (std::numbers::pi / 2.));
+  return Eigen::Matrix4cd{{mc, m0, m0, msm},  // row 0
+                          {m0, m1, m0, m0},   // row 1
+                          {m0, m0, m1, m0},   // row 2
+                          {msp, m0, m0, mc}}; // row 3
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XXPlusYYOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XXPlusYYOp.cpp
index 82f378ab1..07620baba 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XXPlusYYOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/XXPlusYYOp.cpp
@@ -11,13 +11,17 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/Utils/Utils.h"
 
+#include <Eigen/Core>
 #include <cmath>
+#include <complex>
 #include <mlir/Dialect/Arith/IR/Arith.h>
 #include <mlir/IR/Builders.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
 #include <mlir/Support/LogicalResult.h>
+#include <numbers>
+#include <optional>
 #include <variant>
 
 using namespace mlir;
@@ -75,8 +79,9 @@ void XXPlusYYOp::build(OpBuilder& odsBuilder, OperationState& odsState,
                        Value qubit0In, Value qubit1In,
                        const std::variant<double, Value>& theta,
                        const std::variant<double, Value>& beta) {
-  auto thetaOperand = variantToValue(odsBuilder, odsState.location, theta);
-  auto betaOperand = variantToValue(odsBuilder, odsState.location, beta);
+  const auto thetaOperand =
+      variantToValue(odsBuilder, odsState.location, theta);
+  const auto betaOperand = variantToValue(odsBuilder, odsState.location, beta);
   build(odsBuilder, odsState, qubit0In, qubit1In, thetaOperand, betaOperand);
 }
 
@@ -84,3 +89,24 @@ void XXPlusYYOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                              MLIRContext* context) {
   results.add<MergeSubsequentXXPlusYY>(context);
 }
+
+std::optional<Eigen::Matrix4cd> XXPlusYYOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  const auto theta = valueToDouble(getTheta());
+  const auto beta = valueToDouble(getBeta());
+  if (!theta || !beta) {
+    return std::nullopt;
+  }
+
+  const auto m0 = 0.0 + 0i;
+  const auto m1 = 1.0 + 0i;
+  const auto mc = std::cos(*theta / 2.0) + 0i;
+  const auto s = std::sin(*theta / 2.0);
+  const auto msp = std::polar(s, *beta - (std::numbers::pi / 2));
+  const auto msm = std::polar(s, -*beta - (std::numbers::pi / 2));
+  return Eigen::Matrix4cd{{m1, m0, m0, m0},  // row 0
+                          {m0, mc, msm, m0}, // row 1
+                          {m0, msp, mc, m0}, // row 2
+                          {m0, m0, m0, m1}}; // row 3
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/YOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/YOp.cpp
index 4caa205f0..a04f1e4ad 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/YOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/YOp.cpp
@@ -11,6 +11,8 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
+#include <complex>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -39,3 +41,9 @@ void YOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<RemoveSubsequentY>(context);
 }
+
+Eigen::Matrix2cd YOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  return Eigen::Matrix2cd{{0, -1i}, {1i, 0}};
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ZOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ZOp.cpp
index 3ca4334c9..c049449d8 100644
--- a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ZOp.cpp
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/ZOp.cpp
@@ -11,6 +11,7 @@
 #include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/QCOUtils.h"
 
+#include <Eigen/Core>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/IR/OperationSupport.h>
 #include <mlir/IR/PatternMatch.h>
@@ -39,3 +40,8 @@ void ZOp::getCanonicalizationPatterns(RewritePatternSet& results,
                                       MLIRContext* context) {
   results.add<RemoveSubsequentZ>(context);
 }
+
+Eigen::Matrix2cd ZOp::getUnitaryMatrix() {
+  return Eigen::Matrix2cd{{1.0, 0.0},   // row 0
+                          {0.0, -1.0}}; // row 1
+}
diff --git a/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/iSWAPOp.cpp b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/iSWAPOp.cpp
new file mode 100644
index 000000000..846ce142d
--- /dev/null
+++ b/mlir/lib/Dialect/QCO/IR/Operations/StandardGates/iSWAPOp.cpp
@@ -0,0 +1,26 @@
+/*
+ * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
+ * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
+ * All rights reserved.
+ *
+ * SPDX-License-Identifier: MIT
+ *
+ * Licensed under the MIT License
+ */
+
+#include "mlir/Dialect/QCO/IR/QCODialect.h"
+
+#include <Eigen/Core>
+#include <complex>
+
+using namespace mlir;
+using namespace mlir::qco;
+
+Eigen::Matrix4cd iSWAPOp::getUnitaryMatrix() {
+  using namespace std::complex_literals;
+
+  return Eigen::Matrix4cd{{1, 0, 0, 0},  // row 0
+                          {0, 0, 1i, 0}, // row 1
+                          {0, 1i, 0, 0}, // row 2
+                          {0, 0, 0, 1}}; // row 3
+}
diff --git a/mlir/unittests/CMakeLists.txt b/mlir/unittests/CMakeLists.txt
index ed224ed99..86b0e5916 100644
--- a/mlir/unittests/CMakeLists.txt
+++ b/mlir/unittests/CMakeLists.txt
@@ -11,5 +11,6 @@ add_subdirectory(Dialect)
 
 add_custom_target(mqt-core-mlir-unittests)
 
-add_dependencies(mqt-core-mlir-unittests mqt-core-mlir-compiler-pipeline-test
-                 mqt-core-mlir-dialect-qco-ir-modifiers-test mqt-core-mlir-dialect-utils-test)
+add_dependencies(
+  mqt-core-mlir-unittests mqt-core-mlir-compiler-pipeline-test mqt-core-mlir-qco-dialect-test
+  mqt-core-mlir-dialect-qco-ir-modifiers-test mqt-core-mlir-dialect-utils-test)
diff --git a/mlir/unittests/Dialect/QCO/IR/CMakeLists.txt b/mlir/unittests/Dialect/QCO/IR/CMakeLists.txt
index 0a87c6bc8..6e15d1105 100644
--- a/mlir/unittests/Dialect/QCO/IR/CMakeLists.txt
+++ b/mlir/unittests/Dialect/QCO/IR/CMakeLists.txt
@@ -6,4 +6,11 @@
 #
 # Licensed under the MIT License
 
+add_executable(mqt-core-mlir-qco-dialect-test test_unitary_matrix.cpp test_unitary_op_interface.cpp)
+
+target_link_libraries(mqt-core-mlir-qco-dialect-test PRIVATE GTest::gtest_main MLIRQCODialect
+                                                             MLIRLLVMDialect MLIRQCOProgramBuilder)
+
+gtest_discover_tests(mqt-core-mlir-qco-dialect-test)
+
 add_subdirectory(Modifiers)
diff --git a/mlir/unittests/Dialect/QCO/IR/test_unitary_matrix.cpp b/mlir/unittests/Dialect/QCO/IR/test_unitary_matrix.cpp
new file mode 100644
index 000000000..3aa1736a9
--- /dev/null
+++ b/mlir/unittests/Dialect/QCO/IR/test_unitary_matrix.cpp
@@ -0,0 +1,206 @@
+/*
+ * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
+ * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
+ * All rights reserved.
+ *
+ * SPDX-License-Identifier: MIT
+ *
+ * Licensed under the MIT License
+ */
+
+#include "mlir/Dialect/QCO/Builder/QCOProgramBuilder.h"
+#include "mlir/Dialect/QCO/IR/QCODialect.h"
+
+#include <Eigen/Core>
+#include <complex>
+#include <functional>
+#include <gtest/gtest.h>
+#include <llvm/Support/Casting.h>
+#include <llvm/Support/raw_ostream.h>
+#include <memory>
+#include <mlir/Dialect/Arith/IR/Arith.h>
+#include <mlir/Dialect/ControlFlow/IR/ControlFlow.h>
+#include <mlir/Dialect/Func/IR/FuncOps.h>
+#include <mlir/Dialect/LLVMIR/LLVMDialect.h>
+#include <mlir/Dialect/SCF/IR/SCF.h>
+#include <mlir/IR/BuiltinOps.h>
+#include <mlir/IR/DialectRegistry.h>
+#include <mlir/IR/MLIRContext.h>
+#include <mlir/IR/OwningOpRef.h>
+#include <string>
+
+namespace {
+
+using namespace mlir;
+using namespace std::complex_literals;
+
+class QcoUnitaryMatrixTest : public testing::Test {
+protected:
+  void SetUp() override {
+    DialectRegistry registry;
+    registry
+        .insert<qco::QCODialect, arith::ArithDialect, cf::ControlFlowDialect,
+                func::FuncDialect, scf::SCFDialect, LLVM::LLVMDialect>();
+
+    context = std::make_unique<MLIRContext>();
+    context->appendDialectRegistry(registry);
+    context->loadAllAvailableDialects();
+  }
+
+  /**
+   * @brief Build expected QCO IR programmatically and run canonicalization
+   */
+  [[nodiscard]] OwningOpRef<ModuleOp> buildQCOIR(
+      const std::function<void(qco::QCOProgramBuilder&)>& buildFunc) const {
+    qco::QCOProgramBuilder builder(context.get());
+    builder.initialize();
+    buildFunc(builder);
+    return builder.finalize();
+  }
+
+  /**
+   * @brief Get first operation of given type in a module containing a function
+   *        as its first operation.
+   */
+  template <typename OpType>
+  [[nodiscard]] OpType getFirstOp(ModuleOp moduleOp) {
+    auto funcOp = llvm::dyn_cast<func::FuncOp>(
+        moduleOp.getBody()->getOperations().front());
+    if (!funcOp) {
+      return nullptr;
+    }
+
+    auto ops = funcOp.getOps<OpType>();
+    if (ops.empty()) {
+      return nullptr;
+    }
+
+    return *ops.begin();
+  }
+
+  /**
+   * @brief Get text representation of given module.
+   */
+  [[nodiscard]] static std::string toString(ModuleOp moduleOp) {
+    std::string buffer;
+    llvm::raw_string_ostream serializeStream{buffer};
+    moduleOp->print(serializeStream);
+    return serializeStream.str();
+  }
+
+private:
+  std::unique_ptr<MLIRContext> context;
+};
+
+} // namespace
+
+// ##################################################
+// # Standard Gates Unitary Matrix Tests
+// ##################################################
+
+/**
+ * @brief Test: Identity unitary matrix
+ *
+ * @details
+ * Ensure the correct gate definition is returned.
+ */
+TEST_F(QcoUnitaryMatrixTest, IdOpMatrix) {
+  auto moduleOp = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(1, "q");
+    builder.id(reg[0]);
+  });
+  auto op = getFirstOp<qco::IdOp>(*moduleOp);
+  ASSERT_TRUE(op) << toString(*moduleOp);
+
+  EXPECT_EQ(op.getUnitaryMatrix(), Eigen::Matrix2cd::Identity());
+}
+
+/**
+ * @brief Test: X unitary matrix
+ *
+ * @details
+ * Ensure the correct gate definition is returned.
+ */
+TEST_F(QcoUnitaryMatrixTest, XOpMatrix) {
+  auto moduleOp = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(1, "q");
+    builder.x(reg[0]);
+  });
+  auto op = getFirstOp<qco::XOp>(*moduleOp);
+  ASSERT_TRUE(op) << toString(*moduleOp);
+
+  const Eigen::Matrix2cd expectedValue{{0, 1}, {1, 0}};
+  EXPECT_EQ(op.getUnitaryMatrix(), expectedValue);
+}
+
+/**
+ * @brief Test: U2 unitary matrix
+ *
+ * @details
+ * Ensure the correct gate definition is returned.
+ */
+TEST_F(QcoUnitaryMatrixTest, U2OpMatrix) {
+  auto moduleOp = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(1, "q");
+    builder.u2(0.2, 0.8, reg[0]);
+  });
+  auto op = getFirstOp<qco::U2Op>(*moduleOp);
+  ASSERT_TRUE(op) << toString(*moduleOp);
+
+  const Eigen::Matrix2cd expectedValue{
+      {0.70710678 + 0.i, -0.49264604 - 0.50724736i},
+      {0.69301172 + 0.14048043i, // NOLINT(modernize-use-std-numbers)
+       0.38205142 + 0.59500984i}};
+  const auto actualValue = op.getUnitaryMatrix();
+  ASSERT_TRUE(actualValue);
+  EXPECT_TRUE(actualValue->isApprox(expectedValue, 1e-8));
+}
+
+/**
+ * @brief Test: CX unitary matrix
+ *
+ * @details
+ * Ensure the correct gate definition is returned.
+ */
+TEST_F(QcoUnitaryMatrixTest, CtrlXOpMatrix) {
+  auto moduleOp = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(2, "q");
+    builder.cx(reg[0], reg[1]);
+  });
+  auto op = getFirstOp<qco::CtrlOp>(*moduleOp);
+  ASSERT_TRUE(op) << toString(*moduleOp);
+
+  const Eigen::MatrixXcd cxMatrix = Eigen::Matrix4cd{
+      {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}}};
+
+  EXPECT_EQ(op.getUnitaryMatrix(), cxMatrix);
+}
+
+/**
+ * @brief Test: CX unitary matrix, both orientations
+ *
+ * @details
+ * Ensure the correct gate definition is returned and is equal for both
+ * orientations.
+ */
+TEST_F(QcoUnitaryMatrixTest, CtrlX10OpMatrix) {
+  auto moduleOp01 = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(2, "q");
+    builder.cx(reg[0], reg[1]);
+  });
+  auto moduleOp10 = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(2, "q");
+    builder.cx(reg[1], reg[0]);
+  });
+  auto op01 = getFirstOp<qco::CtrlOp>(*moduleOp01);
+  ASSERT_TRUE(op01) << toString(*moduleOp01);
+  auto op10 = getFirstOp<qco::CtrlOp>(*moduleOp10);
+  ASSERT_TRUE(op10) << toString(*moduleOp10);
+
+  const Eigen::MatrixXcd cxMatrix = Eigen::Matrix4cd{
+      {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}}};
+
+  EXPECT_EQ(op01.getUnitaryMatrix(), cxMatrix);
+  EXPECT_EQ(op10.getUnitaryMatrix(), cxMatrix);
+  EXPECT_EQ(op10.getUnitaryMatrix(), op01.getUnitaryMatrix());
+}
diff --git a/mlir/unittests/Dialect/QCO/IR/test_unitary_op_interface.cpp b/mlir/unittests/Dialect/QCO/IR/test_unitary_op_interface.cpp
new file mode 100644
index 000000000..2bfd5cc69
--- /dev/null
+++ b/mlir/unittests/Dialect/QCO/IR/test_unitary_op_interface.cpp
@@ -0,0 +1,320 @@
+/*
+ * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
+ * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
+ * All rights reserved.
+ *
+ * SPDX-License-Identifier: MIT
+ *
+ * Licensed under the MIT License
+ */
+
+#include "mlir/Dialect/QCO/Builder/QCOProgramBuilder.h"
+#include "mlir/Dialect/QCO/IR/QCODialect.h"
+
+#include <array>
+#include <complex>
+#include <cstddef>
+#include <functional>
+#include <gtest/gtest.h>
+#include <llvm/ADT/SmallVector.h>
+#include <llvm/Support/Casting.h>
+#include <llvm/Support/Error.h>
+#include <llvm/Support/raw_ostream.h>
+#include <memory>
+#include <mlir/Dialect/Arith/IR/Arith.h>
+#include <mlir/Dialect/ControlFlow/IR/ControlFlow.h>
+#include <mlir/Dialect/Func/IR/FuncOps.h>
+#include <mlir/Dialect/LLVMIR/LLVMDialect.h>
+#include <mlir/Dialect/SCF/IR/SCF.h>
+#include <mlir/IR/BuiltinOps.h>
+#include <mlir/IR/DialectRegistry.h>
+#include <mlir/IR/MLIRContext.h>
+#include <mlir/IR/OwningOpRef.h>
+#include <string>
+#include <tuple>
+
+namespace {
+
+using namespace mlir;
+using namespace std::complex_literals;
+
+class QcoUnitaryOpInterfaceTest : public testing::Test {
+protected:
+  void SetUp() override {
+    DialectRegistry registry;
+    registry
+        .insert<qco::QCODialect, arith::ArithDialect, cf::ControlFlowDialect,
+                func::FuncDialect, scf::SCFDialect, LLVM::LLVMDialect>();
+
+    context = std::make_unique<MLIRContext>();
+    context->appendDialectRegistry(registry);
+    context->loadAllAvailableDialects();
+  }
+
+  /**
+   * @brief Build expected QCO IR programmatically and run canonicalization
+   */
+  [[nodiscard]] OwningOpRef<ModuleOp> buildQCOIR(
+      const std::function<void(qco::QCOProgramBuilder&)>& buildFunc) const {
+    qco::QCOProgramBuilder builder(context.get());
+    builder.initialize();
+    buildFunc(builder);
+    return builder.finalize();
+  }
+
+  /**
+   * @brief Return list of unitary matrices in given function
+   */
+  template <typename MatrixType>
+  [[nodiscard]] llvm::SmallVector<llvm::Expected<MatrixType>, 0>
+  getMatricesInFunction(func::FuncOp funcOp) {
+    llvm::SmallVector<llvm::Expected<MatrixType>, 0> matrices;
+    for (auto&& op : funcOp.getOps()) {
+      auto unitaryOp = llvm::dyn_cast_if_present<qco::UnitaryOpInterface>(op);
+      if (unitaryOp) {
+        if (auto matrix = unitaryOp.getUnitaryMatrix<MatrixType>()) {
+          matrices.push_back(*matrix);
+        } else {
+          matrices.push_back(llvm::createStringError(
+              "Failed to get matrix of gate '%d' (%s)", matrices.size(),
+              unitaryOp.getBaseSymbol().data()));
+        }
+      }
+    }
+    return matrices;
+  }
+
+  /**
+   * @brief Get text representation of given module.
+   */
+  [[nodiscard]] static std::string toString(ModuleOp moduleOp) {
+    std::string buffer;
+    llvm::raw_string_ostream serializeStream{buffer};
+    moduleOp->print(serializeStream);
+    return serializeStream.str();
+  }
+
+private:
+  std::unique_ptr<MLIRContext> context;
+};
+
+} // namespace
+
+TEST_F(QcoUnitaryOpInterfaceTest, getUnitaryMatrix1x1) {
+  using Matrix1cd = Eigen::Matrix<std::complex<double>, 1, 1>;
+  const auto expectedValue = Matrix1cd{{0.5486898605 + 0.8360259786i}};
+  auto moduleOp =
+      buildQCOIR([](qco::QCOProgramBuilder& builder) { builder.gphase(0.99); });
+
+  auto&& moduleOps = moduleOp->getBody()->getOperations();
+  ASSERT_FALSE(moduleOps.empty());
+  auto funcOp = llvm::dyn_cast<func::FuncOp>(*moduleOps.begin());
+
+  Matrix1cd actualValue = Matrix1cd::Zero();
+  for (auto&& op : funcOp.getOps()) {
+    auto unitaryOp = llvm::dyn_cast<qco::UnitaryOpInterface>(op);
+    if (unitaryOp) {
+      auto matrix = unitaryOp.getUnitaryMatrix<Matrix1cd>();
+      ASSERT_TRUE(matrix) << toString(*moduleOp)
+                          << "\nFailed to get matrix of gate";
+      actualValue = *matrix;
+    }
+  }
+
+  EXPECT_TRUE(actualValue.isApprox(expectedValue, 1e-8))
+      << actualValue << "\nExpected:\n"
+      << expectedValue << "\nDifference:\n"
+      << actualValue - expectedValue;
+}
+
+TEST_F(QcoUnitaryOpInterfaceTest, getUnitaryMatrix2x2) {
+  const auto expectedValues = std::array{
+      Eigen::Matrix2cd{{1, 0}, {0, 1}},
+      Eigen::Matrix2cd{{0.87758256, -0.47942554i}, {-0.47942554i, 0.87758256}},
+      Eigen::Matrix2cd{{0.99500417, -0.09195267 - 0.03887696i},
+                       {0.09537451 + 0.02950279i, 0.76102116 + 0.64099928i}}};
+  auto moduleOp = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(1, "q");
+    reg[0] = builder.id(reg[0]);
+    reg[0] = builder.rx(1.0, reg[0]);
+    reg[0] = builder.u(0.2, 0.3, 0.4, reg[0]);
+  });
+
+  auto&& moduleOps = moduleOp->getBody()->getOperations();
+  ASSERT_FALSE(moduleOps.empty());
+  auto funcOp = llvm::dyn_cast<func::FuncOp>(*moduleOps.begin());
+
+  llvm::SmallVector<Eigen::Matrix2cd> actualValues;
+  for (auto&& op : funcOp.getOps()) {
+    auto unitaryOp = llvm::dyn_cast<qco::UnitaryOpInterface>(op);
+    if (unitaryOp) {
+      auto matrix = unitaryOp.getUnitaryMatrix<Eigen::Matrix2cd>();
+      ASSERT_TRUE(matrix) << toString(*moduleOp)
+                          << "\nFailed to get matrix of gate "
+                          << actualValues.size();
+      actualValues.push_back(*matrix);
+    }
+  }
+
+  ASSERT_EQ(actualValues.size(), expectedValues.size());
+  for (std::size_t i = 0; i < actualValues.size(); ++i) {
+    EXPECT_TRUE(actualValues[i].isApprox(expectedValues.at(i), 1e-8))
+        << "Wrong matrix at gate " << i;
+  }
+}
+
+TEST_F(QcoUnitaryOpInterfaceTest, combine2x2UnitaryMatrices) {
+  // use Qiskit to build same circuit (`qc`) and get unitary using
+  // `qiskit.quantum_info.Operator(qc).data`
+  const auto expectedValue =
+      Eigen::Matrix2cd{{-0.57126014 - 0.1499036i, -0.10275332 - 0.80039522i},
+                       {0.6908436 - 0.41704421i, -0.47252452 - 0.35430187i}};
+  auto moduleOp = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(1, "q");
+    reg[0] = builder.id(reg[0]);
+    reg[0] = builder.h(reg[0]);
+    reg[0] = builder.rx(1.2, reg[0]);
+    reg[0] = builder.ry(1.3, reg[0]);
+    reg[0] = builder.rz(1.4, reg[0]);
+    reg[0] = builder.u(0.5, 0.4, 0.3, reg[0]);
+    reg[0] = builder.x(reg[0]);
+    reg[0] = builder.y(reg[0]);
+    reg[0] = builder.z(reg[0]);
+    reg[0] = builder.sx(reg[0]);
+    reg[0] = builder.s(reg[0]);
+    reg[0] = builder.p(0.2, reg[0]);
+    reg[0] = builder.sdg(reg[0]);
+    reg[0] = builder.t(reg[0]);
+    reg[0] = builder.r(1.0, 0.9, reg[0]);
+    reg[0] = builder.tdg(reg[0]);
+    reg[0] = builder.sxdg(reg[0]);
+
+    // not supported by Qiskit
+    // reg[0] = builder.cgphase(2.7, reg[0]);
+    // reg[0] = builder.u2(2.2, 0.5, reg[0]);
+  });
+
+  auto&& moduleOps = moduleOp->getBody()->getOperations();
+  ASSERT_FALSE(moduleOps.empty());
+  auto funcOp = llvm::dyn_cast<func::FuncOp>(*moduleOps.begin());
+  ASSERT_TRUE(funcOp);
+
+  auto matrices = getMatricesInFunction<Eigen::Matrix2cd>(funcOp);
+
+  Eigen::Matrix2cd combinedMatrix = Eigen::Matrix2cd::Identity();
+  for (auto&& matrix : matrices) {
+    ASSERT_TRUE(static_cast<bool>(matrix))
+        << llvm::toString(matrix.takeError());
+    combinedMatrix = *matrix * combinedMatrix;
+  }
+
+  EXPECT_TRUE(combinedMatrix.isApprox(expectedValue, 1e-8))
+      << "Combination of matrices does not match expected matrix";
+}
+
+TEST_F(QcoUnitaryOpInterfaceTest, combine4x4UnitaryMatrices) {
+  // use Qiskit to build same circuit (`qc`) and get unitary using
+  // `qiskit.quantum_info.Operator(qc).data`:
+  // qc = QuantumCircuit(2, 0)
+  // qc.rxx(1.1, 0, 1)
+  // qc.ryy(1.2, 0, 1)
+  // qc.swap(0, 1)
+  // qc.rzz(1.3, 0, 1)
+  // qc.rzx(1.4, 1, 0)
+  // qc.iswap(0, 1)
+  // qc.ecr(1, 0)
+  // qc.dcx(1, 0)
+  // qc.append(XXMinusYYGate(2.0, 2.1), [0, 1])
+  // qc.append(XXPlusYYGate(2.2, 2.3), [0, 1])
+  // qc.cx(1, 0)
+  // qc.ch(1, 0)
+  // print(Operator(qc).data)
+
+  const auto expectedValue =
+      Eigen::Matrix4cd{{-0.19081581 - 0.2947213i, 0.20121632 - 0.46723087i,
+                        -0.37097846 + 0.25410653i, 0.01790367 + 0.64453107i},
+                       {0.67932253 - 0.04649638i, -0.11526485 - 0.61915801i,
+                        -0.0866555 - 0.22376126i, 0.20493411 - 0.20034995i},
+                       {0.57996486 - 0.00202388i, -0.23429021 + 0.41790602i,
+                        0.20814033 + 0.16778316i, 0.06442165 + 0.5987283i},
+                       {0.14845576 + 0.23384873i, 0.05966105 + 0.33386807i,
+                        -0.78994584 + 0.21151491i, 0.33101926 - 0.17297858i}};
+  auto moduleOp = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(2, "q");
+
+    std::tie(reg[0], reg[1]) = builder.rxx(1.1, reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.ryy(1.2, reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.swap(reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.rzz(1.3, reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.rzx(1.4, reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.iswap(reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.ecr(reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.dcx(reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.xx_minus_yy(2.0, 2.1, reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.xx_plus_yy(2.2, 2.3, reg[0], reg[1]);
+
+    // implicit conversion of dynamic matrix to fixed-size matrix,
+    // if size matches; note: Qiskit respects qubit order, but QCO does not have
+    // concept of an "order" between qubits
+    std::tie(reg[0], reg[1]) = builder.cx(reg[1], reg[0]);
+    std::tie(reg[0], reg[1]) = builder.ch(reg[1], reg[0]);
+  });
+
+  auto&& moduleOps = moduleOp->getBody()->getOperations();
+  ASSERT_FALSE(moduleOps.empty());
+  auto funcOp = llvm::dyn_cast<func::FuncOp>(*moduleOps.begin());
+  ASSERT_TRUE(funcOp);
+
+  auto matrices = getMatricesInFunction<Eigen::Matrix4cd>(funcOp);
+
+  Eigen::Matrix4cd combinedMatrix = Eigen::Matrix4cd::Identity();
+  for (auto&& matrix : matrices) {
+    ASSERT_TRUE(static_cast<bool>(matrix))
+        << llvm::toString(matrix.takeError());
+    combinedMatrix = *matrix * combinedMatrix;
+  }
+
+  EXPECT_TRUE(combinedMatrix.isApprox(expectedValue, 1e-8))
+      << "Combination of matrices does not match expected matrix\nCombined:\n"
+      << combinedMatrix << "\nExpected:\n"
+      << expectedValue << "\nDifference:\n"
+      << combinedMatrix - expectedValue;
+}
+
+TEST_F(QcoUnitaryOpInterfaceTest, getDynamicUnitaryMatrix) {
+  const auto expectedValues = std::array{
+      Eigen::MatrixXcd{{1, 0}, {0, 1}},
+      Eigen::MatrixXcd{{1, 0, 0, 0}, {0, 0, 1, 0}, {0, 1, 0, 0}, {0, 0, 0, 1}},
+      Eigen::MatrixXcd{{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}},
+      Eigen::MatrixXcd{{1, 0, 0, 0, 0, 0, 0, 0},
+                       {0, 1, 0, 0, 0, 0, 0, 0},
+                       {0, 0, 1, 0, 0, 0, 0, 0},
+                       {0, 0, 0, 1, 0, 0, 0, 0},
+                       {0, 0, 0, 0, 1, 0, 0, 0},
+                       {0, 0, 0, 0, 0, 1, 0, 0},
+                       {0, 0, 0, 0, 0, 0, 0, 1},
+                       {0, 0, 0, 0, 0, 0, 1, 0}},
+  };
+  auto moduleOp = buildQCOIR([](qco::QCOProgramBuilder& builder) {
+    auto reg = builder.allocQubitRegister(3, "q");
+    reg[0] = builder.id(reg[0]);
+    std::tie(reg[0], reg[1]) = builder.swap(reg[0], reg[1]);
+    std::tie(reg[0], reg[1]) = builder.cx(reg[0], reg[1]);
+    std::ignore = builder.mcx(std::array{reg[0], reg[1]}, reg[2]);
+  });
+
+  auto&& moduleOps = moduleOp->getBody()->getOperations();
+  ASSERT_FALSE(moduleOps.empty());
+  auto funcOp = llvm::dyn_cast<func::FuncOp>(*moduleOps.begin());
+  ASSERT_TRUE(funcOp);
+
+  auto actualValues = getMatricesInFunction<Eigen::MatrixXcd>(funcOp);
+
+  ASSERT_EQ(actualValues.size(), expectedValues.size());
+  for (std::size_t i = 0; i < actualValues.size(); ++i) {
+    ASSERT_TRUE(static_cast<bool>(actualValues[i]))
+        << llvm::toString(actualValues[i].takeError());
+    EXPECT_TRUE(actualValues[i]->isApprox(expectedValues.at(i), 1e-8))
+        << "Wrong matrix at gate " << i;
+  }
+}