diff --git a/Project.toml b/Project.toml
index 6521764f..c49b76c2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "LatticeTools"
 uuid = "8a267902-5e0e-44c3-8ef8-9b2b5c352816"
 author = "Kyungmin Lee"
-version = "0.8.1"
+version = "0.9.0"
 
 [deps]
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
diff --git a/src/Irrep/irrep.jl b/src/Irrep/irrep.jl
index 5e1a5f0b..734ed6f1 100644
--- a/src/Irrep/irrep.jl
+++ b/src/Irrep/irrep.jl
@@ -7,7 +7,7 @@ export group_order
 export get_irrep_components
 export get_irrep_iterator
 export little_symmetry
-
+export make_product_irrep
 
 struct IrrepData
     group::FiniteGroup
@@ -212,6 +212,23 @@ function iscompatible(
 end
 
 
+
+function make_product_irrep(irrep_iterators...; tol::Real=Base.rtoldefault(Float64))
+    nested_symops_and_amplitude_list = [
+        [(x, y) for (x, y) in irrep_iterator if !isapprox(y, zero(y); atol=tol)]
+            for irrep_iterator in irrep_iterators
+    ]
+    symops_and_amplitudes = [
+        (prod(reverse([op for (op, amp) in elems])), prod(amp for (op, amp) in elems))
+            for elems in Iterators.product(nested_symops_and_amplitude_list...)
+            # elems has the form of ((s, phis), (t, phit), (p, phit), ...)
+            # we want the resulting element to be (p*t*s, ϕp*ϕt*ϕs). the phase commutes, while operations do not necessarily. 
+    ]
+    return symops_and_amplitudes
+end
+
+
+
 # struct SymmorphicSpaceIrrepComponent{S1<:SymmetryOrEmbedding, S2<:SymmetryOrEmbedding} <:AbstractSymmetryIrrepComponent
 
 #     component1::IrrepComponent{S1} # e.g. Translation
diff --git a/src/LatticeTools.jl b/src/LatticeTools.jl
index e1c002fe..5abb202d 100644
--- a/src/LatticeTools.jl
+++ b/src/LatticeTools.jl
@@ -14,6 +14,7 @@ include("SymmetryOperation/identityoperation.jl")
 include("SymmetryOperation/translationoperation.jl")
 include("SymmetryOperation/pointoperation.jl")
 include("SymmetryOperation/spaceoperation.jl")
+include("SymmetryOperation/productoperation.jl")
 
 include("Symmetry/abstractsymmetry.jl")
 include("Symmetry/translationsymmetry.jl")
diff --git a/src/SymmetryEmbedding/sitepermutation.jl b/src/SymmetryEmbedding/sitepermutation.jl
index c9d58952..7467d645 100644
--- a/src/SymmetryEmbedding/sitepermutation.jl
+++ b/src/SymmetryEmbedding/sitepermutation.jl
@@ -7,7 +7,7 @@ export embed
 export isidentity
 
 
-abstract type AbstractSymmetryOperationEmbedding end
+abstract type AbstractSymmetryOperationEmbedding <: AbstractSymmetryOperation end
 abstract type AbstractSpaceSymmetryOperationEmbedding <: AbstractSymmetryOperationEmbedding end
 
 
@@ -19,7 +19,7 @@ Represents a permutation of sites as a symmetry operation of a lattice.
 # Fields
 * `permutation::Permutation`
 """
-struct SitePermutation <:AbstractSpaceSymmetryOperationEmbedding
+struct SitePermutation <: AbstractSpaceSymmetryOperationEmbedding
     permutation::Permutation
     SitePermutation(p::Permutation) = new(p)
     SitePermutation(p::AbstractVector{<:Integer}) = new(Permutation(p))
diff --git a/src/SymmetryOperation/productoperation.jl b/src/SymmetryOperation/productoperation.jl
index 64344d27..b913bd6d 100644
--- a/src/SymmetryOperation/productoperation.jl
+++ b/src/SymmetryOperation/productoperation.jl
@@ -1,3 +1,54 @@
+export DirectProductOperation
+export isidentity
+import LinearAlgebra.×
+export directproduct
+
+"""
+    DirectProductOperation{Ops}
+
+Direct product of symmetry operations.
+math```
+  g = (g_1, g_2, \\ldots)
+```
+"""
+struct DirectProductOperation{Ops<:Tuple{Vararg{AbstractSymmetryOperation}}}<:AbstractSymmetryOperation
+    operations::Ops
+    function DirectProductOperation(ops::Vararg{AbstractSymmetryOperation})
+        T = typeof(ops)
+        return new{T}(ops)
+    end
+    function DirectProductOperation(ops::T) where {T<:Tuple{Vararg{AbstractSymmetryOperation}}}
+        return new{T}(ops)
+    end
+end
+
+×(lhs::AbstractSymmetryOperation, rhs::AbstractSymmetryOperation) = DirectProductOperation(lhs, rhs)
+×(lhs::DirectProductOperation, rhs::AbstractSymmetryOperation) = DirectProductOperation(lhs.operations..., rhs)
+×(lhs::AbstractSymmetryOperation, rhs::DirectProductOperation) = DirectProductOperation(lhs, rhs.operations...)
+×(lhs::DirectProductOperation, rhs::DirectProductOperation) = DirectProductOperation(lhs.operations..., rhs.operations...)
+
+Base.hash(x::P, h::UInt) where {P<:DirectProductOperation} = Base.hash(P, Base.hash(x.operations, h))
+
+Base.:(*)(lhs::P, rhs::P) where {P<:DirectProductOperation} = DirectProductOperation(lhs.operations .* rhs.operations)
+Base.:(^)(obj::DirectProductOperation, n::Integer) = DirectProductOperation(obj.operations.^n)
+Base.:(==)(lhs::P, rhs::P) where {P<:DirectProductOperation} = all(lhs.operations .== rhs.operations)
+
+Base.inv(obj::DirectProductOperation) = DirectProductOperation(Base.inv.(obj.operations))
+
+function Base.isapprox(lhs::P, rhs::P; atol::Real=0, rtol::Real=Base.rtoldefault(Float64)) where {P<:DirectProductOperation}
+    return all(isapprox.(lhs.operations, rhs.operations; atol=atol, rtol=rtol))
+end
+
+isidentity(obj::DirectProductOperation) = all(isidentity, obj.operations)
+
+function directproduct(::Type{E}, products::Function...) where {E<:DirectProductOperation}
+    function product(lhs::E, rhs::E)
+        return DirectProductOperation([p(l, r) for (p, l, r) in zip(products, lhs.operations, rhs.operations)]...)
+    end
+    return product
+end
+
+
 # export ProductOperation
 
 # export apply_operation