diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0b426e88b..3ce866819 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Rename `sparse_to_dense` as `to_dense` and `dense_to_sparse` as `to_sparse`. ([#392])
 - Fix erroneous definition of the stochastic term in `smesolve`. ([#393])
+- Change name of `MultiSiteOperator` to `multisite_operator`. ([#394])
 
 ## [v0.26.0]
 Release date: 2025-02-09
@@ -120,3 +121,4 @@ Release date: 2024-11-13
 [#389]: https://github.com/qutip/QuantumToolbox.jl/issues/389
 [#392]: https://github.com/qutip/QuantumToolbox.jl/issues/392
 [#393]: https://github.com/qutip/QuantumToolbox.jl/issues/393
+[#394]: https://github.com/qutip/QuantumToolbox.jl/issues/394
diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md
index 1a5af52a8..6e95ec396 100644
--- a/docs/src/resources/api.md
+++ b/docs/src/resources/api.md
@@ -260,7 +260,7 @@ fidelity
 
 ```@docs
 Lattice
-MultiSiteOperator
+multisite_operator
 DissipativeIsing
 ```
 
diff --git a/docs/src/users_guide/cluster.md b/docs/src/users_guide/cluster.md
index 3e1726ffa..f290ae1b4 100644
--- a/docs/src/users_guide/cluster.md
+++ b/docs/src/users_guide/cluster.md
@@ -99,9 +99,9 @@ hy = 0.0
 hz = 0.0
 γ = 1
 
-Sx = mapreduce(i -> MultiSiteOperator(latt, i=>sigmax()), +, 1:latt.N)
-Sy = mapreduce(i -> MultiSiteOperator(latt, i=>sigmay()), +, 1:latt.N)
-Sz = mapreduce(i -> MultiSiteOperator(latt, i=>sigmaz()), +, 1:latt.N)
+Sx = mapreduce(i -> multisite_operator(latt, i=>sigmax()), +, 1:latt.N)
+Sy = mapreduce(i -> multisite_operator(latt, i=>sigmay()), +, 1:latt.N)
+Sz = mapreduce(i -> multisite_operator(latt, i=>sigmaz()), +, 1:latt.N)
 
 H, c_ops = DissipativeIsing(Jx, Jy, Jz, hx, hy, hz, γ, latt; boundary_condition = Val(:periodic_bc), order = 1)
 e_ops = [Sx, Sy, Sz]
diff --git a/src/deprecated.jl b/src/deprecated.jl
index a5d729711..56e5d04e7 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -93,3 +93,7 @@ correlation_2op_1t(
 } = error(
     "The parameter order of `correlation_2op_1t` has been changed, please use `?correlation_2op_1t` to check the updated docstring.",
 )
+
+MultiSiteOperator(dims::Union{AbstractVector,Tuple}, pairs::Pair{<:Integer,<:QuantumObject}...) = error(
+    "`MultiSiteOperator` has been deprecated and will be removed in next major release, please use `multisite_operator` instead.",
+)
diff --git a/src/spin_lattice.jl b/src/spin_lattice.jl
index 99f356df3..17275aba1 100644
--- a/src/spin_lattice.jl
+++ b/src/spin_lattice.jl
@@ -1,4 +1,4 @@
-export Lattice, MultiSiteOperator, DissipativeIsing
+export Lattice, multisite_operator, DissipativeIsing
 
 @doc raw"""
     Lattice
@@ -15,7 +15,7 @@ end
 
 #Definition of many-body operators
 @doc raw"""
-    MultiSiteOperator(dims::Union{AbstractVector, Tuple}, pairs::Pair{<:Integer,<:QuantumObject}...)
+    multisite_operator(dims::Union{AbstractVector, Tuple}, pairs::Pair{<:Integer,<:QuantumObject}...)
 
 A Julia function for generating a multi-site operator ``\\hat{O} = \\hat{O}_i \\hat{O}_j \\cdots \\hat{O}_k``. The function takes a vector of dimensions `dims` and a list of pairs `pairs` where the first element of the pair is the site index and the second element is the operator acting on that site.
 
@@ -28,7 +28,7 @@ A Julia function for generating a multi-site operator ``\\hat{O} = \\hat{O}_i \\
 
 # Example
 ```jldoctest
-julia> op = MultiSiteOperator(Val(8), 5=>sigmax(), 7=>sigmaz());
+julia> op = multisite_operator(Val(8), 5=>sigmax(), 7=>sigmaz());
 
 julia> op.dims
 8-element SVector{8, Int64} with indices SOneTo(8):
@@ -42,7 +42,7 @@ julia> op.dims
  2
 ```
 """
-function MultiSiteOperator(dims::Union{AbstractVector,Tuple}, pairs::Pair{<:Integer,<:QuantumObject}...)
+function multisite_operator(dims::Union{AbstractVector,Tuple}, pairs::Pair{<:Integer,<:QuantumObject}...)
     sites_unsorted = collect(first.(pairs))
     idxs = sortperm(sites_unsorted)
     _sites = sites_unsorted[idxs]
@@ -61,13 +61,13 @@ function MultiSiteOperator(dims::Union{AbstractVector,Tuple}, pairs::Pair{<:Inte
 
     return QuantumObject(data; type = Operator, dims = dims)
 end
-function MultiSiteOperator(N::Union{Integer,Val}, pairs::Pair{<:Integer,<:QuantumObject}...)
+function multisite_operator(N::Union{Integer,Val}, pairs::Pair{<:Integer,<:QuantumObject}...)
     dims = ntuple(j -> 2, makeVal(N))
 
-    return MultiSiteOperator(dims, pairs...)
+    return multisite_operator(dims, pairs...)
 end
-function MultiSiteOperator(latt::Lattice, pairs::Pair{<:Integer,<:QuantumObject}...)
-    return MultiSiteOperator(makeVal(latt.N), pairs...)
+function multisite_operator(latt::Lattice, pairs::Pair{<:Integer,<:QuantumObject}...)
+    return multisite_operator(makeVal(latt.N), pairs...)
 end
 
 #Definition of nearest-neighbour sites on lattice
@@ -127,7 +127,7 @@ function DissipativeIsing(
     boundary_condition::Union{Symbol,Val} = Val(:periodic_bc),
     order::Integer = 1,
 )
-    S = [MultiSiteOperator(latt, i => sigmam()) for i in 1:latt.N]
+    S = [multisite_operator(latt, i => sigmam()) for i in 1:latt.N]
     c_ops = sqrt(γ) .* S
 
     op_sum(S, i::CartesianIndex) =
@@ -135,17 +135,17 @@ function DissipativeIsing(
 
     H = 0
     if (Jx != 0 || hx != 0)
-        S = [MultiSiteOperator(latt, i => sigmax()) for i in 1:latt.N]
+        S = [multisite_operator(latt, i => sigmax()) for i in 1:latt.N]
         H += Jx / 2 * mapreduce(i -> op_sum(S, i), +, latt.car_idx) #/2 because we are double counting
         H += hx * sum(S)
     end
     if (Jy != 0 || hy != 0)
-        S = [MultiSiteOperator(latt, i => sigmay()) for i in 1:latt.N]
+        S = [multisite_operator(latt, i => sigmay()) for i in 1:latt.N]
         H += Jy / 2 * mapreduce(i -> op_sum(S, i), +, latt.car_idx)
         H += hy * sum(S)
     end
     if (Jz != 0 || hz != 0)
-        S = [MultiSiteOperator(latt, i => sigmaz()) for i in 1:latt.N]
+        S = [multisite_operator(latt, i => sigmaz()) for i in 1:latt.N]
         H += Jz / 2 * mapreduce(i -> op_sum(S, i), +, latt.car_idx)
         H += hz * sum(S)
     end
diff --git a/test/core-test/low_rank_dynamics.jl b/test/core-test/low_rank_dynamics.jl
index c7ce03c06..7b718c734 100644
--- a/test/core-test/low_rank_dynamics.jl
+++ b/test/core-test/low_rank_dynamics.jl
@@ -15,12 +15,12 @@
     i = 1
     for j in 1:N_modes
         i += 1
-        i <= M && (ϕ[i] = MultiSiteOperator(latt, j => sigmap()) * ϕ[1])
+        i <= M && (ϕ[i] = multisite_operator(latt, j => sigmap()) * ϕ[1])
     end
     for k in 1:N_modes-1
         for l in k+1:N_modes
             i += 1
-            i <= M && (ϕ[i] = MultiSiteOperator(latt, k => sigmap(), l => sigmap()) * ϕ[1])
+            i <= M && (ϕ[i] = multisite_operator(latt, k => sigmap(), l => sigmap()) * ϕ[1])
         end
     end
     for i in i+1:M
@@ -43,11 +43,11 @@
     hz = 0.0
     γ = 1
 
-    Sx = mapreduce(i -> MultiSiteOperator(latt, i => sigmax()), +, 1:latt.N)
-    Sy = mapreduce(i -> MultiSiteOperator(latt, i => sigmay()), +, 1:latt.N)
-    Sz = mapreduce(i -> MultiSiteOperator(latt, i => sigmaz()), +, 1:latt.N)
+    Sx = mapreduce(i -> multisite_operator(latt, i => sigmax()), +, 1:latt.N)
+    Sy = mapreduce(i -> multisite_operator(latt, i => sigmay()), +, 1:latt.N)
+    Sz = mapreduce(i -> multisite_operator(latt, i => sigmaz()), +, 1:latt.N)
     SFxx = mapreduce(
-        x -> MultiSiteOperator(latt, x[1] => sigmax(), x[2] => sigmax()),
+        x -> multisite_operator(latt, x[1] => sigmax(), x[2] => sigmax()),
         +,
         Iterators.product(1:latt.N, 1:latt.N),
     )