diff --git a/qiskit/quantum_info/states/quantum_state.py b/qiskit/quantum_info/states/quantum_state.py
index 2e2bcb672649..640f7fec3c99 100644
--- a/qiskit/quantum_info/states/quantum_state.py
+++ b/qiskit/quantum_info/states/quantum_state.py
@@ -454,7 +454,7 @@ def _subsystem_probabilities(probs, dims, qargs=None):
         if qargs is None:
             return probs
         # Convert qargs to tensor axes
-        probs_tens = np.reshape(probs, dims)
+        probs_tens = np.reshape(probs, list(reversed(dims)))
         ndim = probs_tens.ndim
         qargs_axes = [ndim - 1 - i for i in reversed(qargs)]
         # Get sum axis for marginalized subsystems
diff --git a/releasenotes/notes/probabilities_dict_bug_fix-aac3b3d3853828dc.yaml b/releasenotes/notes/probabilities_dict_bug_fix-aac3b3d3853828dc.yaml
new file mode 100644
index 000000000000..e6c2590cfee7
--- /dev/null
+++ b/releasenotes/notes/probabilities_dict_bug_fix-aac3b3d3853828dc.yaml
@@ -0,0 +1,7 @@
+---
+fixes:
+  - |
+    Fixed an issue with the :meth:`.Statevector.probabilities_dict` and :meth:`.DensityMatrix.probabilities_dict`
+    method where it would return incorrect results for non-qubit systems when the ``qargs`` argument was
+    specified.
+    Fixed `#9210 <https://github.com/Qiskit/qiskit-terra/issues/9210>`__
\ No newline at end of file
diff --git a/test/python/quantum_info/states/test_densitymatrix.py b/test/python/quantum_info/states/test_densitymatrix.py
index aad19dc9e2da..6a619094ce3e 100644
--- a/test/python/quantum_info/states/test_densitymatrix.py
+++ b/test/python/quantum_info/states/test_densitymatrix.py
@@ -678,6 +678,67 @@ def test_sample_counts_w(self):
                 counts = state.sample_counts(shots, qargs=qargs)
                 self.assertDictAlmostEqual(counts, target, threshold)
 
+    def test_probabilities_dict_unequal_dims(self):
+        """Test probabilities_dict for a state with unequal subsystem dimensions."""
+
+        vec = np.zeros(60, dtype=float)
+        vec[15:20] = np.ones(5)
+        vec[40:46] = np.ones(6)
+        state = DensityMatrix(vec / np.sqrt(11.0), dims=[3, 4, 5])
+
+        p = 1.0 / 11.0
+
+        self.assertDictEqual(
+            state.probabilities_dict(),
+            {
+                s: p
+                for s in [
+                    "110",
+                    "111",
+                    "112",
+                    "120",
+                    "121",
+                    "311",
+                    "312",
+                    "320",
+                    "321",
+                    "322",
+                    "330",
+                ]
+            },
+        )
+
+        # differences due to rounding
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[0]), {"0": 4 * p, "1": 4 * p, "2": 3 * p}, delta=1e-10
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[1]), {"1": 5 * p, "2": 5 * p, "3": p}, delta=1e-10
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[2]), {"1": 5 * p, "3": 6 * p}, delta=1e-10
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[0, 1]),
+            {"10": p, "11": 2 * p, "12": 2 * p, "20": 2 * p, "21": 2 * p, "22": p, "30": p},
+            delta=1e-10,
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[1, 0]),
+            {"01": p, "11": 2 * p, "21": 2 * p, "02": 2 * p, "12": 2 * p, "22": p, "03": p},
+            delta=1e-10,
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[0, 2]),
+            {"10": 2 * p, "11": 2 * p, "12": p, "31": 2 * p, "32": 2 * p, "30": 2 * p},
+            delta=1e-10,
+        )
+
     def test_sample_counts_qutrit(self):
         """Test sample_counts method for qutrit state"""
         p = 0.3
diff --git a/test/python/quantum_info/states/test_statevector.py b/test/python/quantum_info/states/test_statevector.py
index 876bd1c6684c..976e31372e79 100644
--- a/test/python/quantum_info/states/test_statevector.py
+++ b/test/python/quantum_info/states/test_statevector.py
@@ -710,6 +710,67 @@ def test_sample_counts_w(self):
                 counts = state.sample_counts(shots, qargs=qargs)
                 self.assertDictAlmostEqual(counts, target, threshold)
 
+    def test_probabilities_dict_unequal_dims(self):
+        """Test probabilities_dict for a state with unequal subsystem dimensions."""
+
+        vec = np.zeros(60, dtype=float)
+        vec[15:20] = np.ones(5)
+        vec[40:46] = np.ones(6)
+        state = Statevector(vec / np.sqrt(11.0), dims=[3, 4, 5])
+
+        p = 1.0 / 11.0
+
+        self.assertDictEqual(
+            state.probabilities_dict(),
+            {
+                s: p
+                for s in [
+                    "110",
+                    "111",
+                    "112",
+                    "120",
+                    "121",
+                    "311",
+                    "312",
+                    "320",
+                    "321",
+                    "322",
+                    "330",
+                ]
+            },
+        )
+
+        # differences due to rounding
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[0]), {"0": 4 * p, "1": 4 * p, "2": 3 * p}, delta=1e-10
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[1]), {"1": 5 * p, "2": 5 * p, "3": p}, delta=1e-10
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[2]), {"1": 5 * p, "3": 6 * p}, delta=1e-10
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[0, 1]),
+            {"10": p, "11": 2 * p, "12": 2 * p, "20": 2 * p, "21": 2 * p, "22": p, "30": p},
+            delta=1e-10,
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[1, 0]),
+            {"01": p, "11": 2 * p, "21": 2 * p, "02": 2 * p, "12": 2 * p, "22": p, "03": p},
+            delta=1e-10,
+        )
+
+        self.assertDictAlmostEqual(
+            state.probabilities_dict(qargs=[0, 2]),
+            {"10": 2 * p, "11": 2 * p, "12": p, "31": 2 * p, "32": 2 * p, "30": 2 * p},
+            delta=1e-10,
+        )
+
     def test_sample_counts_qutrit(self):
         """Test sample_counts method for qutrit state"""
         p = 0.3