diff --git a/qiskit/algorithms/eigen_solvers/eigen_solver.py b/qiskit/algorithms/eigen_solvers/eigen_solver.py
index 5f34d10e0c16..fb46137459bb 100644
--- a/qiskit/algorithms/eigen_solvers/eigen_solver.py
+++ b/qiskit/algorithms/eigen_solvers/eigen_solver.py
@@ -13,12 +13,16 @@
 """The Eigensolver interface"""
 
 from abc import ABC, abstractmethod
-from typing import List, Optional
+from typing import Dict, Optional, List, Union, TypeVar
 
 import numpy as np
 from qiskit.opflow import OperatorBase
 from ..algorithm_result import AlgorithmResult
 
+# Introduced new type to maintain readability.
+_T = TypeVar("_T")  # Pylint does not allow single character class names.
+ListOrDict = Union[List[Optional[_T]], Dict[str, _T]]
+
 
 class Eigensolver(ABC):
     """The Eigensolver Interface.
@@ -30,7 +34,7 @@ class Eigensolver(ABC):
 
     @abstractmethod
     def compute_eigenvalues(
-        self, operator: OperatorBase, aux_operators: Optional[List[Optional[OperatorBase]]] = None
+        self, operator: OperatorBase, aux_operators: Optional[ListOrDict[OperatorBase]] = None
     ) -> "EigensolverResult":
         """
         Computes eigenvalues. Operator and aux_operators can be supplied here and
@@ -90,11 +94,11 @@ def eigenstates(self, value: np.ndarray) -> None:
         self._eigenstates = value
 
     @property
-    def aux_operator_eigenvalues(self) -> Optional[np.ndarray]:
+    def aux_operator_eigenvalues(self) -> Optional[List[ListOrDict[complex]]]:
         """return aux operator eigen values"""
         return self._aux_operator_eigenvalues
 
     @aux_operator_eigenvalues.setter
-    def aux_operator_eigenvalues(self, value: np.ndarray) -> None:
+    def aux_operator_eigenvalues(self, value: List[ListOrDict[complex]]) -> None:
         """set aux operator eigen values"""
         self._aux_operator_eigenvalues = value
diff --git a/qiskit/algorithms/eigen_solvers/numpy_eigen_solver.py b/qiskit/algorithms/eigen_solvers/numpy_eigen_solver.py
index 3e8434973d65..54a2b81a70a5 100755
--- a/qiskit/algorithms/eigen_solvers/numpy_eigen_solver.py
+++ b/qiskit/algorithms/eigen_solvers/numpy_eigen_solver.py
@@ -20,9 +20,8 @@
 
 from qiskit.opflow import I, ListOp, OperatorBase, StateFn
 from qiskit.utils.validation import validate_min
-
 from ..exceptions import AlgorithmError
-from .eigen_solver import Eigensolver, EigensolverResult
+from .eigen_solver import Eigensolver, EigensolverResult, ListOrDict
 
 logger = logging.getLogger(__name__)
 
@@ -47,7 +46,7 @@ def __init__(
         self,
         k: int = 1,
         filter_criterion: Callable[
-            [Union[List, np.ndarray], float, Optional[List[float]]], bool
+            [Union[List, np.ndarray], float, Optional[ListOrDict[float]]], bool
         ] = None,
     ) -> None:
         """
@@ -85,7 +84,7 @@ def k(self, k: int) -> None:
     @property
     def filter_criterion(
         self,
-    ) -> Optional[Callable[[Union[List, np.ndarray], float, Optional[List[float]]], bool]]:
+    ) -> Optional[Callable[[Union[List, np.ndarray], float, Optional[ListOrDict[float]]], bool]]:
         """returns the filter criterion if set"""
         return self._filter_criterion
 
@@ -93,7 +92,7 @@ def filter_criterion(
     def filter_criterion(
         self,
         filter_criterion: Optional[
-            Callable[[Union[List, np.ndarray], float, Optional[List[float]]], bool]
+            Callable[[Union[List, np.ndarray], float, Optional[ListOrDict[float]]], bool]
         ],
     ) -> None:
         """set the filter criterion"""
@@ -150,7 +149,7 @@ def _get_ground_state_energy(self, operator: OperatorBase) -> None:
             self._solve(operator)
 
     def _get_energies(
-        self, operator: OperatorBase, aux_operators: Optional[List[OperatorBase]]
+        self, operator: OperatorBase, aux_operators: Optional[ListOrDict[OperatorBase]]
     ) -> None:
         if self._ret.eigenvalues is None or self._ret.eigenstates is None:
             self._solve(operator)
@@ -165,12 +164,19 @@ def _get_energies(
 
     @staticmethod
     def _eval_aux_operators(
-        aux_operators: List[OperatorBase], wavefn, threshold: float = 1e-12
-    ) -> np.ndarray:
-        values = []  # type: List[Tuple[float, int]]
-        for operator in aux_operators:
+        aux_operators: ListOrDict[OperatorBase], wavefn, threshold: float = 1e-12
+    ) -> ListOrDict[Tuple[float, int]]:
+
+        # As a list, aux_operators can contain None operators for which None values are returned.
+        # As a dict, the None operators in aux_operators have been dropped in compute_eigenvalues.
+        if isinstance(aux_operators, list):
+            values = [None] * len(aux_operators)  # type: ListOrDict[Tuple[float, int]]
+            key_op_iterator = enumerate(aux_operators)
+        else:
+            values = {}
+            key_op_iterator = aux_operators.items()
+        for key, operator in key_op_iterator:
             if operator is None:
-                values.append(None)
                 continue
             value = 0.0
             if operator.coeff != 0:
@@ -183,11 +189,11 @@ def _eval_aux_operators(
                 else:
                     value = StateFn(operator, is_measurement=True).eval(wavefn)
                 value = value.real if abs(value.real) > threshold else 0.0
-            values.append((value, 0))
-        return np.array(values, dtype=object)
+            values[key] = (value, 0)
+        return values
 
     def compute_eigenvalues(
-        self, operator: OperatorBase, aux_operators: Optional[List[Optional[OperatorBase]]] = None
+        self, operator: OperatorBase, aux_operators: Optional[ListOrDict[OperatorBase]] = None
     ) -> EigensolverResult:
         super().compute_eigenvalues(operator, aux_operators)
 
@@ -195,10 +201,16 @@ def compute_eigenvalues(
             raise AlgorithmError("Operator was never provided")
 
         self._check_set_k(operator)
-        if aux_operators:
-            zero_op = I.tensorpower(operator.num_qubits) * 0.0
+        zero_op = I.tensorpower(operator.num_qubits) * 0.0
+        if isinstance(aux_operators, list) and len(aux_operators) > 0:
             # For some reason Chemistry passes aux_ops with 0 qubits and paulis sometimes.
             aux_operators = [zero_op if op == 0 else op for op in aux_operators]
+        elif isinstance(aux_operators, dict) and len(aux_operators) > 0:
+            aux_operators = {
+                key: zero_op if op == 0 else op  # Convert zero values to zero operators
+                for key, op in aux_operators.items()
+                if op is not None  # Discard None values
+            }
         else:
             aux_operators = None
 
diff --git a/qiskit/algorithms/minimum_eigen_solvers/minimum_eigen_solver.py b/qiskit/algorithms/minimum_eigen_solvers/minimum_eigen_solver.py
index 3c96dbae3ca8..62335b9715c0 100644
--- a/qiskit/algorithms/minimum_eigen_solvers/minimum_eigen_solver.py
+++ b/qiskit/algorithms/minimum_eigen_solvers/minimum_eigen_solver.py
@@ -13,12 +13,16 @@
 """The Minimum Eigensolver interface"""
 
 from abc import ABC, abstractmethod
-from typing import List, Optional
+from typing import Dict, Optional, List, Union, TypeVar
 
 import numpy as np
 from qiskit.opflow import OperatorBase
 from ..algorithm_result import AlgorithmResult
 
+# Introduced new type to maintain readability.
+_T = TypeVar("_T")  # Pylint does not allow single character class names.
+ListOrDict = Union[List[Optional[_T]], Dict[str, _T]]
+
 
 class MinimumEigensolver(ABC):
     """The Minimum Eigensolver Interface.
@@ -30,7 +34,7 @@ class MinimumEigensolver(ABC):
 
     @abstractmethod
     def compute_minimum_eigenvalue(
-        self, operator: OperatorBase, aux_operators: Optional[List[Optional[OperatorBase]]] = None
+        self, operator: OperatorBase, aux_operators: Optional[ListOrDict[OperatorBase]] = None
     ) -> "MinimumEigensolverResult":
         """
         Computes minimum eigenvalue. Operator and aux_operators can be supplied here and
@@ -94,11 +98,11 @@ def eigenstate(self, value: np.ndarray) -> None:
         self._eigenstate = value
 
     @property
-    def aux_operator_eigenvalues(self) -> Optional[np.ndarray]:
+    def aux_operator_eigenvalues(self) -> Optional[ListOrDict[complex]]:
         """return aux operator eigen values"""
         return self._aux_operator_eigenvalues
 
     @aux_operator_eigenvalues.setter
-    def aux_operator_eigenvalues(self, value: np.ndarray) -> None:
+    def aux_operator_eigenvalues(self, value: ListOrDict[complex]) -> None:
         """set aux operator eigen values"""
         self._aux_operator_eigenvalues = value
diff --git a/qiskit/algorithms/minimum_eigen_solvers/numpy_minimum_eigen_solver.py b/qiskit/algorithms/minimum_eigen_solvers/numpy_minimum_eigen_solver.py
index 1d2b29fd12cc..369485976cdf 100644
--- a/qiskit/algorithms/minimum_eigen_solvers/numpy_minimum_eigen_solver.py
+++ b/qiskit/algorithms/minimum_eigen_solvers/numpy_minimum_eigen_solver.py
@@ -18,7 +18,7 @@
 
 from qiskit.opflow import OperatorBase
 from ..eigen_solvers.numpy_eigen_solver import NumPyEigensolver
-from .minimum_eigen_solver import MinimumEigensolver, MinimumEigensolverResult
+from .minimum_eigen_solver import MinimumEigensolver, MinimumEigensolverResult, ListOrDict
 
 logger = logging.getLogger(__name__)
 
@@ -31,7 +31,7 @@ class NumPyMinimumEigensolver(MinimumEigensolver):
     def __init__(
         self,
         filter_criterion: Callable[
-            [Union[List, np.ndarray], float, Optional[List[float]]], bool
+            [Union[List, np.ndarray], float, Optional[ListOrDict[float]]], bool
         ] = None,
     ) -> None:
         """
@@ -49,7 +49,7 @@ def __init__(
     @property
     def filter_criterion(
         self,
-    ) -> Optional[Callable[[Union[List, np.ndarray], float, Optional[List[float]]], bool]]:
+    ) -> Optional[Callable[[Union[List, np.ndarray], float, Optional[ListOrDict[float]]], bool]]:
         """returns the filter criterion if set"""
         return self._ces.filter_criterion
 
@@ -57,7 +57,7 @@ def filter_criterion(
     def filter_criterion(
         self,
         filter_criterion: Optional[
-            Callable[[Union[List, np.ndarray], float, Optional[List[float]]], bool]
+            Callable[[Union[List, np.ndarray], float, Optional[ListOrDict[float]]], bool]
         ],
     ) -> None:
         """set the filter criterion"""
@@ -68,7 +68,7 @@ def supports_aux_operators(cls) -> bool:
         return NumPyEigensolver.supports_aux_operators()
 
     def compute_minimum_eigenvalue(
-        self, operator: OperatorBase, aux_operators: Optional[List[Optional[OperatorBase]]] = None
+        self, operator: OperatorBase, aux_operators: Optional[ListOrDict[OperatorBase]] = None
     ) -> MinimumEigensolverResult:
         super().compute_minimum_eigenvalue(operator, aux_operators)
         result_ces = self._ces.compute_eigenvalues(operator, aux_operators)
diff --git a/qiskit/algorithms/minimum_eigen_solvers/vqe.py b/qiskit/algorithms/minimum_eigen_solvers/vqe.py
index 32224abff65c..9da7012f2ec6 100755
--- a/qiskit/algorithms/minimum_eigen_solvers/vqe.py
+++ b/qiskit/algorithms/minimum_eigen_solvers/vqe.py
@@ -42,7 +42,7 @@
 from qiskit.utils import QuantumInstance, algorithm_globals
 from ..optimizers import Optimizer, SLSQP
 from ..variational_algorithm import VariationalAlgorithm, VariationalResult
-from .minimum_eigen_solver import MinimumEigensolver, MinimumEigensolverResult
+from .minimum_eigen_solver import MinimumEigensolver, MinimumEigensolverResult, ListOrDict
 from ..exceptions import AlgorithmError
 
 logger = logging.getLogger(__name__)
@@ -411,32 +411,42 @@ def supports_aux_operators(cls) -> bool:
     def _eval_aux_ops(
         self,
         parameters: np.ndarray,
-        aux_operators: List[OperatorBase],
+        aux_operators: ListOrDict[OperatorBase],
         expectation: ExpectationBase,
         threshold: float = 1e-12,
-    ) -> np.ndarray:
+    ) -> ListOrDict[complex]:
         # Create new CircuitSampler to avoid breaking existing one's caches.
         sampler = CircuitSampler(self.quantum_instance)
 
-        aux_op_meas = expectation.convert(StateFn(ListOp(aux_operators), is_measurement=True))
+        if isinstance(aux_operators, dict):
+            list_op = ListOp(list(aux_operators.values()))
+        else:
+            list_op = ListOp(aux_operators)
+
+        aux_op_meas = expectation.convert(StateFn(list_op, is_measurement=True))
         aux_op_expect = aux_op_meas.compose(CircuitStateFn(self.ansatz.bind_parameters(parameters)))
         values = np.real(sampler.convert(aux_op_expect).eval())
 
         # Discard values below threshold
         aux_op_results = values * (np.abs(values) > threshold)
-        # Deal with the aux_op behavior where there can be Nones or Zero qubit Paulis in the list
-        _aux_op_nones = [op is None for op in aux_operators]
-        aux_operator_eigenvalues = [
-            None if is_none else [result]
-            for (is_none, result) in zip(_aux_op_nones, aux_op_results)
-        ]
-        # As this has mixed types, since it can included None, it needs to explicitly pass object
-        # data type to avoid numpy 1.19 warning message about implicit conversion being deprecated
-        aux_operator_eigenvalues = np.array([aux_operator_eigenvalues], dtype=object)
+
+        # Return None eigenvalues for None operators if aux_operators is a list.
+        # None operators are already dropped in compute_minimum_eigenvalue if aux_operators is a dict.
+        if isinstance(aux_operators, list):
+            aux_operator_eigenvalues = [None] * len(aux_operators)
+            key_value_iterator = enumerate(aux_op_results)
+        else:
+            aux_operator_eigenvalues = {}
+            key_value_iterator = zip(aux_operators.keys(), aux_op_results)
+
+        for key, value in key_value_iterator:
+            if aux_operators[key] is not None:
+                aux_operator_eigenvalues[key] = value
+
         return aux_operator_eigenvalues
 
     def compute_minimum_eigenvalue(
-        self, operator: OperatorBase, aux_operators: Optional[List[Optional[OperatorBase]]] = None
+        self, operator: OperatorBase, aux_operators: Optional[ListOrDict[OperatorBase]] = None
     ) -> MinimumEigensolverResult:
         super().compute_minimum_eigenvalue(operator, aux_operators)
 
@@ -454,19 +464,24 @@ def compute_minimum_eigenvalue(
         initial_point = _validate_initial_point(self.initial_point, self.ansatz)
 
         bounds = _validate_bounds(self.ansatz)
-
-        # We need to handle the array entries being Optional i.e. having value None
+        # We need to handle the array entries being zero or Optional i.e. having value None
         if aux_operators:
             zero_op = I.tensorpower(operator.num_qubits) * 0.0
-            converted = []
-            for op in aux_operators:
-                if op is None:
-                    converted.append(zero_op)
-                else:
-                    converted.append(op)
 
-            # For some reason Chemistry passes aux_ops with 0 qubits and paulis sometimes.
-            aux_operators = [zero_op if op == 0 else op for op in converted]
+            # Convert the None and zero values when aux_operators is a list.
+            # Drop None and convert zero values when aux_operators is a dict.
+            if isinstance(aux_operators, list):
+                key_op_iterator = enumerate(aux_operators)
+                converted = [zero_op] * len(aux_operators)
+            else:
+                key_op_iterator = aux_operators.items()
+                converted = {}
+            for key, op in key_op_iterator:
+                if op is not None:
+                    converted[key] = zero_op if op == 0 else op
+
+            aux_operators = converted
+
         else:
             aux_operators = None
 
@@ -532,7 +547,7 @@ def compute_minimum_eigenvalue(
 
         if aux_operators is not None:
             aux_values = self._eval_aux_ops(opt_result.x, aux_operators, expectation=expectation)
-            result.aux_operator_eigenvalues = aux_values[0]
+            result.aux_operator_eigenvalues = aux_values
 
         return result
 
diff --git a/releasenotes/notes/support-dict-for-aux-operators-c3c9ad380c208afd.yaml b/releasenotes/notes/support-dict-for-aux-operators-c3c9ad380c208afd.yaml
new file mode 100644
index 000000000000..a487e76601a0
--- /dev/null
+++ b/releasenotes/notes/support-dict-for-aux-operators-c3c9ad380c208afd.yaml
@@ -0,0 +1,9 @@
+---
+features:
+  - |
+    The ``EigenSolver`` and ``MinimumEigenSolver`` interfaces now support the type
+    ``Dict[str, Optional[OperatorBase]]`` for the ``aux_operators`` parameter in the respective
+    ``compute_eigenvalues`` and ``compute_minimum_eigenvalue`` methods.
+    In this case, the auxiliary eigenvalues are also stored in a dictionary under the same keys
+    provided by the `aux_operators` dictionary. Keys that correspond to an operator that does not commute
+    with the main operator are dropped.
\ No newline at end of file
diff --git a/test/python/algorithms/test_numpy_minimum_eigen_solver.py b/test/python/algorithms/test_numpy_minimum_eigen_solver.py
index 98842ef12cef..14ce9f86b09c 100644
--- a/test/python/algorithms/test_numpy_minimum_eigen_solver.py
+++ b/test/python/algorithms/test_numpy_minimum_eigen_solver.py
@@ -40,14 +40,16 @@ def setUp(self):
 
         aux_op1 = PauliSumOp.from_list([("II", 2.0)])
         aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
-        self.aux_ops = [aux_op1, aux_op2]
+        self.aux_ops_list = [aux_op1, aux_op2]
+        self.aux_ops_dict = {"aux_op1": aux_op1, "aux_op2": aux_op2}
 
     def test_cme(self):
         """Basic test"""
         algo = NumPyMinimumEigensolver()
-        result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
-        self.assertEqual(result.eigenvalue.dtype, np.float64)
-        self.assertAlmostEqual(result.eigenvalue, -1.85727503)
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_list
+        )
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[0], [2, 0])
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[1], [0, 0])
@@ -62,9 +64,10 @@ def test_cme_reuse(self):
         self.assertIsNone(result.aux_operator_eigenvalues)
 
         # Add aux_operators and go again
-        result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
-        self.assertEqual(result.eigenvalue.dtype, np.float64)
-        self.assertAlmostEqual(result.eigenvalue, -1.85727503)
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_list
+        )
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[0], [2, 0])
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[1], [0, 0])
@@ -76,15 +79,16 @@ def test_cme_reuse(self):
         self.assertIsNone(result.aux_operator_eigenvalues)
 
         # Set aux_operators and go again
-        result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
-        self.assertEqual(result.eigenvalue.dtype, np.float64)
-        self.assertAlmostEqual(result.eigenvalue, -1.85727503)
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_list
+        )
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[0], [2, 0])
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[1], [0, 0])
 
         # Finally just set one of aux_operators and main operator, remove aux_operators
-        result = algo.compute_minimum_eigenvalue(operator=self.aux_ops[0], aux_operators=[])
+        result = algo.compute_minimum_eigenvalue(operator=self.aux_ops_list[0], aux_operators=[])
         self.assertAlmostEqual(result.eigenvalue, 2 + 0j)
         self.assertIsNone(result.aux_operator_eigenvalues)
 
@@ -97,9 +101,10 @@ def criterion(x, v, a_v):
             return v >= -0.5
 
         algo = NumPyMinimumEigensolver(filter_criterion=criterion)
-        result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
-        self.assertEqual(result.eigenvalue.dtype, np.float64)
-        self.assertAlmostEqual(result.eigenvalue, -0.22491125)
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_list
+        )
+        self.assertAlmostEqual(result.eigenvalue, -0.22491125 + 0j)
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[0], [2, 0])
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[1], [0, 0])
@@ -113,7 +118,9 @@ def criterion(x, v, a_v):
             return False
 
         algo = NumPyMinimumEigensolver(filter_criterion=criterion)
-        result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_list
+        )
         self.assertEqual(result.eigenvalue, None)
         self.assertEqual(result.eigenstate, None)
         self.assertEqual(result.aux_operator_eigenvalues, None)
@@ -125,6 +132,32 @@ def test_cme_1q(self, op):
         result = algo.compute_minimum_eigenvalue(operator=op)
         self.assertAlmostEqual(result.eigenvalue, -1)
 
+    def test_cme_aux_ops_dict(self):
+        """Test dictionary compatibility of aux_operators"""
+        # Start with an empty dictionary
+        algo = NumPyMinimumEigensolver()
+        result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators={})
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertIsNone(result.aux_operator_eigenvalues)
+
+        # Add aux_operators dictionary and go again
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_dict
+        )
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertEqual(len(result.aux_operator_eigenvalues), 2)
+        np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues["aux_op1"], [2, 0])
+        np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues["aux_op2"], [0, 0])
+
+        # Add None and zero operators and go again
+        extra_ops = {"None_op": None, "zero_op": 0, **self.aux_ops_dict}
+        result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=extra_ops)
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertEqual(len(result.aux_operator_eigenvalues), 3)
+        np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues["aux_op1"], [2, 0])
+        np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues["aux_op2"], [0, 0])
+        self.assertEqual(result.aux_operator_eigenvalues["zero_op"], (0.0, 0))
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/python/algorithms/test_vqe.py b/test/python/algorithms/test_vqe.py
index 3a84acfec7c8..aecc7189be24 100644
--- a/test/python/algorithms/test_vqe.py
+++ b/test/python/algorithms/test_vqe.py
@@ -473,6 +473,35 @@ def wrapped_run(circuits, **kwargs):
 
         self.assertEqual(callcount["count"], expected)
 
+    def test_aux_operators_dict(self):
+        """Test dictionary compatibility of aux_operators"""
+        wavefunction = self.ry_wavefunction
+        vqe = VQE(ansatz=wavefunction, quantum_instance=self.statevector_simulator)
+
+        # Start with an empty dictionary
+        result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators={})
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=6)
+        self.assertIsNone(result.aux_operator_eigenvalues)
+
+        # Go again with two auxiliary operators
+        aux_op1 = PauliSumOp.from_list([("II", 2.0)])
+        aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+        aux_ops = {"aux_op1": aux_op1, "aux_op2": aux_op2}
+        result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators=aux_ops)
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=6)
+        self.assertEqual(len(result.aux_operator_eigenvalues), 2)
+        self.assertAlmostEqual(result.aux_operator_eigenvalues["aux_op1"], 2, places=6)
+        self.assertAlmostEqual(result.aux_operator_eigenvalues["aux_op2"], 0, places=6)
+
+        # Go again with additional None and zero operators
+        extra_ops = {**aux_ops, "None_operator": None, "zero_operator": 0}
+        result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators=extra_ops)
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=6)
+        self.assertEqual(len(result.aux_operator_eigenvalues), 3)
+        self.assertAlmostEqual(result.aux_operator_eigenvalues["aux_op1"], 2, places=6)
+        self.assertAlmostEqual(result.aux_operator_eigenvalues["aux_op2"], 0, places=6)
+        self.assertEqual(result.aux_operator_eigenvalues["zero_operator"], 0.0)
+
 
 if __name__ == "__main__":
     unittest.main()