diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0985357d3..b22a252ab 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Changes to `SteadyStateODESolver`. ([#537])
   - Introduce the tolerances for `steadystate` terminate condition (two new fields: `terminate_reltol = 1e-5` and `terminate_abstol = 1e-7`)
   - Fix keyword argument handling for `SteadyStateODESolver` before passing to `mesolve`.
+- Fix incorrect `negativity` and `partial_transpose` for arbitrary subsystem dimension. ([#539])
 
 ## [v0.34.1]
 Release date: 2025-08-23
@@ -308,3 +309,4 @@ Release date: 2024-11-13
 [#531]: https://github.com/qutip/QuantumToolbox.jl/issues/531
 [#536]: https://github.com/qutip/QuantumToolbox.jl/issues/536
 [#537]: https://github.com/qutip/QuantumToolbox.jl/issues/537
+[#539]: https://github.com/qutip/QuantumToolbox.jl/issues/539
diff --git a/src/negativity.jl b/src/negativity.jl
index 77e5d37d9..286ce37a0 100644
--- a/src/negativity.jl
+++ b/src/negativity.jl
@@ -85,20 +85,20 @@ end
 
 # for dense matrices
 function _partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
-    mask2 = [1 + Int(i) for i in mask]
+    nsys = length(mask)
+    mask2 = reverse([mask[s] ? 2 : 1 for s in 1:nsys])
     # mask2 has elements with values equal to 1 or 2
-    #   1 - the subsystem don't need to be transposed
-    #   2 - the subsystem need be transposed
+    #   1 - the subsystem (in reversed order) don't need to be transposed
+    #   2 - the subsystem (in reversed order) need to be transposed
 
-    nsys = length(mask2)
-    dims = dimensions_to_dims(get_dimensions_to(ρ))
+    dims_rev = reverse(dimensions_to_dims(get_dimensions_to(ρ)))
     pt_dims = reshape(Vector(1:(2*nsys)), (nsys, 2))
     pt_idx = [
         [pt_dims[n, mask2[n]] for n in 1:nsys]   # origin   value in mask2
         [pt_dims[n, 3-mask2[n]] for n in 1:nsys] # opposite value in mask2 (1 -> 2, and 2 -> 1)
     ]
     return QuantumObject(
-        reshape(permutedims(reshape(ρ.data, (dims..., dims...)), pt_idx), size(ρ)),
+        reshape(permutedims(reshape(ρ.data, (dims_rev..., dims_rev...)), pt_idx), size(ρ)),
         Operator(),
         Dimensions(ρ.dimensions.to),
     )
@@ -110,7 +110,8 @@ function _partial_transpose(
     mask::Vector{Bool},
 ) where {DimsType<:AbstractDimensions}
     M, N = size(ρ)
-    dimsTuple = Tuple(dimensions_to_dims(get_dimensions_to(ρ)))
+    dims_rev = reverse(Tuple(dimensions_to_dims(get_dimensions_to(ρ))))
+    mask_rev = reverse(mask)
     colptr = ρ.data.colptr
     rowval = ρ.data.rowval
     nzval = ρ.data.nzval
@@ -130,13 +131,13 @@ function _partial_transpose(
                 I_pt[n] = i
                 J_pt[n] = j
             else
-                ket_pt = [Base._ind2sub(dimsTuple, i)...]
-                bra_pt = [Base._ind2sub(dimsTuple, j)...]
-                for sys in findall(m -> m, mask)
+                ket_pt = [Base._ind2sub(dims_rev, i)...]
+                bra_pt = [Base._ind2sub(dims_rev, j)...]
+                for sys in findall(m -> m, mask_rev)
                     @inbounds ket_pt[sys], bra_pt[sys] = bra_pt[sys], ket_pt[sys]
                 end
-                I_pt[n] = Base._sub2ind(dimsTuple, ket_pt...)
-                J_pt[n] = Base._sub2ind(dimsTuple, bra_pt...)
+                I_pt[n] = Base._sub2ind(dims_rev, ket_pt...)
+                J_pt[n] = Base._sub2ind(dims_rev, bra_pt...)
             end
             V_pt[n] = nzval[p]
         end
diff --git a/test/core-test/negativity_and_partial_transpose.jl b/test/core-test/negativity_and_partial_transpose.jl
index 44ad9deb6..a554d6065 100644
--- a/test/core-test/negativity_and_partial_transpose.jl
+++ b/test/core-test/negativity_and_partial_transpose.jl
@@ -1,6 +1,6 @@
 @testitem "Negativity and Partial Transpose" begin
     @testset "negativity" begin
-        rho = (1 / 40) * Qobj(
+        rho1 = (1 / 40) * Qobj(
             [
                 15 1 1 15
                 1 5 -3 1
@@ -9,15 +9,32 @@
             ];
             dims = (2, 2),
         )
-        Neg = negativity(rho, 1)
-        @test Neg ≈ 0.25
-        @test negativity(rho, 2) ≈ Neg
-        @test negativity(rho, 1; logarithmic = true) ≈ log2(2 * Neg + 1)
-        @test_throws ArgumentError negativity(rho, 3)
+        Neg1 = negativity(rho1, 1)
+        @test Neg1 ≈ 0.25
+        @test negativity(rho1, 2) ≈ Neg1
+        @test negativity(rho1, 1; logarithmic = true) ≈ log2(2 * Neg1 + 1)
+        @test_throws ArgumentError negativity(rho1, 3)
+
+        # a maximally entanglment state with subsystem dimension (3, 2):
+        # (|1,0⟩ - i|2,1⟩) / √2
+        rho2_d = ket2dm((tensor(basis(3, 1), basis(2, 0)) - 1im * tensor(basis(3, 2), basis(2, 1))) / √2)
+        rho2_s = to_sparse(rho2_d)
+        @test negativity(rho2_d, 1) ≈ 0.5
+        @test negativity(rho2_d, 2) ≈ 0.5
+        @test negativity(rho2_s, 1) ≈ 0.5
+        @test negativity(rho2_s, 2) ≈ 0.5
+
+        # a separable state with subsystem dimension (3, 2)
+        rho3_d = tensor(rand_dm(3), rand_dm(2))
+        rho3_s = to_sparse(rho3_d)
+        @test abs(negativity(rho3_d, 1)) < 1e-10
+        @test abs(negativity(rho3_d, 2)) < 1e-10
+        @test abs(negativity(rho3_s, 1)) < 1e-10
+        @test abs(negativity(rho3_s, 2)) < 1e-10
 
         @testset "Type Inference (negativity)" begin
-            @inferred negativity(rho, 1)
-            @inferred negativity(rho, 1; logarithmic = true)
+            @inferred negativity(rho1, 1)
+            @inferred negativity(rho1, 1; logarithmic = true)
         end
     end