diff --git a/.pylintdict b/.pylintdict
index 245f309e..6846bce4 100644
--- a/.pylintdict
+++ b/.pylintdict
@@ -362,6 +362,7 @@ trainability
 transpilation
 transpile
 transpiled
+transpiler
 trotterization
 trotterized
 uncompute
diff --git a/TODO.md b/TODO.md
new file mode 100644
index 00000000..a29cf8f3
--- /dev/null
+++ b/TODO.md
@@ -0,0 +1,39 @@
+- [x] phase_estimators/hamiltonian_phase_estimation.py
+- [x] phase_estimators/ipe.py
+- [x] phase_estimators/phase_estimation.py
+- [x] eigensolvers/vqd.py
+- [x] amplitude_amplifiers/grover.py (test doesn't work but should probably be refactored or removed)
+- [x] time_evolvers/pvqd/utils.py
+- [x] time_evolvers/pvqd/pvqd.py
+- [x] time_evolvers/trotterization/trotter_qrte.py
+- [x] time_evolvers/variational/variational_principles/imaginary_mc_lachlan_principle.py
+- [x] time_evolvers/variational/variational_principles/real_mc_lachlan_principle.py
+- [x] time_evolvers/variational/var_qite.py
+- [x] time_evolvers/variational/var_qrte.py
+- [x] time_evolvers/variational/var_qte.py
+- [x] state_fidelities/compute_uncompute.py
+- [x] optimizers/qnspsa.py
+- [x] optimizers/umda.py
+- [x] optimizers/spsa.py
+- [x] observables_evaluator.py
+- [] gradients/reverse/reverse_gradient.py
+- [] gradients/reverse/reverse_qgt.py
+- [] gradients/finite_diff/finite_diff_estimator_gradient.py
+- [] gradients/finite_diff/finite_diff_sampler_gradient.py
+- [] gradients/spsa/spsa_estimator_gradient.py
+- [] gradients/spsa/spsa_sampler_gradient.py
+- [] gradients/lin_comb/lin_comb_sampler_gradient.py
+- [] gradients/lin_comb/lin_comb_estimator_gradient.py
+- [] gradients/lin_comb/lin_comb_qgt.py
+- [] gradients/base/base_sampler_gradient.py
+- [] gradients/base/base_qgt.py
+- [] gradients/base/base_estimator_gradient.py
+- [x] minimum_eigensolvers/vqe.py
+- [x] minimum_eigensolvers/adapt_vqe.py
+- [x] minimum_eigensolvers/qaoa.py
+- [x] minimum_eigensolvers/diagonal_estimator.py
+- [x] minimum_eigensolvers/sampling_vqe.py
+- [] amplitude_estimators/mlae.py
+- [] amplitude_estimators/fae.py
+- [] amplitude_estimators/iae.py
+- [] amplitude_estimators/ae.py
diff --git a/docs/tutorials/10_pvqd.ipynb b/docs/tutorials/10_pvqd.ipynb
index 239276fe..463a3174 100644
--- a/docs/tutorials/10_pvqd.ipynb
+++ b/docs/tutorials/10_pvqd.ipynb
@@ -125,7 +125,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from qiskit.primitives import Sampler, Estimator\n",
+    "from qiskit.primitives import Sampler, StatevectorEstimator as Estimator\n",
     "from qiskit_algorithms.state_fidelities import ComputeUncompute\n",
     "\n",
     "# the fidelity is used to evaluate the objective: the overlap of the variational form and the trotter step\n",
diff --git a/qiskit_algorithms/algorithm_job.py b/qiskit_algorithms/algorithm_job.py
index b1215574..0ec841bc 100644
--- a/qiskit_algorithms/algorithm_job.py
+++ b/qiskit_algorithms/algorithm_job.py
@@ -20,26 +20,6 @@ class AlgorithmJob(PrimitiveJob):
     """
     This class is introduced for typing purposes and provides no
     additional function beyond that inherited from its parents.
-
-    Update: :meth:`AlgorithmJob.submit()` method added. See its
-    documentation for more info.
     """
 
-    def submit(self) -> None:
-        """
-        Submit the job for execution.
-
-        For V1 primitives, Qiskit ``PrimitiveJob`` subclassed JobV1 and defined ``submit()``.
-        ``PrimitiveJob`` was updated for V2 primitives, no longer subclasses ``JobV1``, and
-        now has a private ``_submit()`` method, with ``submit()`` being deprecated as of
-        Qiskit version 0.46. This maintains the ``submit()`` for ``AlgorithmJob`` here as
-        it's called in many places for such a job. An alternative could be to make
-        0.46 the required minimum version and alter all algorithm's call sites to use
-        ``_submit()`` and make this an empty class again as it once was. For now this
-        way maintains compatibility with the current min version of 0.44.
-        """
-        # TODO: Considering changing this in the future - see above docstring.
-        try:
-            super()._submit()
-        except AttributeError:
-            super().submit()
+    pass
diff --git a/qiskit_algorithms/amplitude_amplifiers/grover.py b/qiskit_algorithms/amplitude_amplifiers/grover.py
index 56fd2ad1..aaab79d6 100644
--- a/qiskit_algorithms/amplitude_amplifiers/grover.py
+++ b/qiskit_algorithms/amplitude_amplifiers/grover.py
@@ -18,16 +18,14 @@
 from typing import Any
 
 import numpy as np
-
 from qiskit import ClassicalRegister, QuantumCircuit
-from qiskit.primitives import BaseSampler
-from qiskit.quantum_info import Statevector
+from qiskit.primitives import BaseSamplerV2
 
 from qiskit_algorithms.exceptions import AlgorithmError
 from qiskit_algorithms.utils import algorithm_globals
-
 from .amplification_problem import AmplificationProblem
 from .amplitude_amplifier import AmplitudeAmplifier, AmplitudeAmplifierResult
+from ..custom_types import Transpiler
 
 
 class Grover(AmplitudeAmplifier):
@@ -116,7 +114,9 @@ def __init__(
         iterations: list[int] | Iterator[int] | int | None = None,
         growth_rate: float | None = None,
         sample_from_iterations: bool = False,
-        sampler: BaseSampler | None = None,
+        sampler: BaseSamplerV2 | None = None,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
     ) -> None:
         r"""
         Args:
@@ -136,6 +136,11 @@ def __init__(
                 powers of the Grover operator, a random integer sample between 0 and smaller value
                 than the iteration is used as a power, see [1], Section 4.
             sampler: A Sampler to use for sampling the results of the circuits.
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced within this algorithm. If set to `None`, these won't be
+                transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
 
         Raises:
             ValueError: If ``growth_rate`` is a float but not larger than 1.
@@ -165,9 +170,11 @@ def __init__(
         self._sampler = sampler
         self._sample_from_iterations = sample_from_iterations
         self._iterations_arg = iterations
+        self._transpiler = transpiler
+        self._transpiler_options = transpiler_options if transpiler_options is not None else {}
 
     @property
-    def sampler(self) -> BaseSampler | None:
+    def sampler(self) -> BaseSamplerV2 | None:
         """Get the sampler.
 
         Returns:
@@ -176,7 +183,7 @@ def sampler(self) -> BaseSampler | None:
         return self._sampler
 
     @sampler.setter
-    def sampler(self, sampler: BaseSampler) -> None:
+    def sampler(self, sampler: BaseSamplerV2) -> None:
         """Set the sampler.
 
         Args:
@@ -234,23 +241,29 @@ def amplify(self, amplification_problem: AmplificationProblem) -> "GroverResult"
             # sample from [0, power) if specified
             if self._sample_from_iterations:
                 power = algorithm_globals.random.integers(power)
+
             # Run a grover experiment for a given power of the Grover operator.
-            if self._sampler is not None:
-                qc = self.construct_circuit(amplification_problem, power, measurement=True)
-                job = self._sampler.run([qc])
-
-                try:
-                    results = job.result()
-                except Exception as exc:
-                    raise AlgorithmError("Sampler job failed.") from exc
-
-                num_bits = len(amplification_problem.objective_qubits)
-                circuit_results: dict[str, Any] | Statevector | np.ndarray = {
-                    np.binary_repr(k, num_bits): v for k, v in results.quasi_dists[0].items()
-                }
-                top_measurement, max_probability = max(
-                    circuit_results.items(), key=lambda x: x[1]  # type: ignore[union-attr]
-                )
+            qc = self.construct_circuit(amplification_problem, power, measurement=True)
+
+            if self._transpiler is not None:
+                qc = self._transpiler.run(qc, **self._transpiler_options)
+
+            job = self._sampler.run([qc])
+
+            try:
+                results = job.result()
+            except Exception as exc:
+                raise AlgorithmError("Sampler job failed.") from exc
+
+            circuit_results = getattr(results[0].data, qc.cregs[0].name)
+            circuit_results = {
+                label: value / circuit_results.num_shots
+                for label, value in circuit_results.get_counts().items()
+            }
+
+            top_measurement, max_probability = max(
+                circuit_results.items(), key=lambda x: x[1]  # type: ignore[union-attr]
+            )
 
             all_circuit_results.append(circuit_results)
 
diff --git a/qiskit_algorithms/custom_types.py b/qiskit_algorithms/custom_types.py
new file mode 100644
index 00000000..3d98fd3e
--- /dev/null
+++ b/qiskit_algorithms/custom_types.py
@@ -0,0 +1,28 @@
+# This code is part of a Qiskit project.
+#
+# (C) Copyright IBM 2024.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Types used by the qiskit-algorithms package."""
+from __future__ import annotations
+
+from typing import Any, Protocol, Union
+
+from qiskit import QuantumCircuit
+
+_Circuits = Union[list[QuantumCircuit], QuantumCircuit]
+
+
+class Transpiler(Protocol):
+    """A Generic type to represent a transpiler."""
+
+    def run(self, circuits: _Circuits, **options: Any) -> _Circuits:
+        """Transpile a circuit or a list of quantum circuits."""
+        pass
diff --git a/qiskit_algorithms/eigensolvers/vqd.py b/qiskit_algorithms/eigensolvers/vqd.py
index 48ba25bb..a0d4d420 100644
--- a/qiskit_algorithms/eigensolvers/vqd.py
+++ b/qiskit_algorithms/eigensolvers/vqd.py
@@ -25,7 +25,7 @@
 import numpy as np
 
 from qiskit.circuit import QuantumCircuit
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 from qiskit.quantum_info import SparsePauliOp
 
@@ -88,7 +88,7 @@ class VQD(VariationalAlgorithm, Eigensolver):
     updated once the VQD object has been constructed.
 
     Attributes:
-            estimator (BaseEstimator): The primitive instance used to perform the expectation
+            estimator (BaseEstimatorV2): The primitive instance used to perform the expectation
                 estimation as indicated in the VQD paper.
             fidelity (BaseStateFidelity): The fidelity class instance used to compute the
                 overlap estimation as indicated in the VQD paper.
@@ -115,7 +115,7 @@ class VQD(VariationalAlgorithm, Eigensolver):
 
     def __init__(
         self,
-        estimator: BaseEstimator,
+        estimator: BaseEstimatorV2,
         fidelity: BaseStateFidelity,
         ansatz: QuantumCircuit,
         optimizer: Optimizer | Minimizer | Sequence[Optimizer | Minimizer],
@@ -363,9 +363,9 @@ def compute_eigenvalues(
 
                     raise AlgorithmError(
                         f"Convergence threshold is set to {self.convergence_threshold} but an "
-                        f"average fidelity of {average_fidelity:.5f} with the previous eigenstates"
-                        f"has been observed during the evaluation of the {step}{suffix} lowest"
-                        f"eigenvalue."
+                        f"average (weighted by the betas) fidelity of {average_fidelity:.5f} with "
+                        f"the previous eigenstates has been observed during the evaluation of the "
+                        f"{step}{suffix} lowest eigenvalue."
                     )
                 logger.info(
                     (
@@ -433,9 +433,7 @@ def evaluate_energy(parameters: np.ndarray) -> float | np.ndarray:
                 parameters = np.reshape(parameters, (-1, num_parameters))
             batch_size = len(parameters)
 
-            estimator_job = self.estimator.run(
-                batch_size * [self.ansatz], batch_size * [operator], parameters
-            )
+            estimator_job = self.estimator.run([(self.ansatz, operator, parameters)])
 
             total_cost = np.zeros(batch_size)
 
@@ -454,18 +452,17 @@ def evaluate_energy(parameters: np.ndarray) -> float | np.ndarray:
                     total_cost += np.real(betas[state] * cost)
 
             try:
-                estimator_result = estimator_job.result()
+                estimator_result = estimator_job.result()[0]
 
             except Exception as exc:
                 raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc
 
-            values = estimator_result.values + total_cost
+            values = estimator_result.data.evs + total_cost
 
             if self.callback is not None:
-                metadata = estimator_result.metadata
-                for params, value, meta in zip(parameters, values, metadata):
+                for params, value in zip(parameters, values):
                     self._eval_count += 1
-                    self.callback(self._eval_count, params, value, meta, step)
+                    self.callback(self._eval_count, params, value, estimator_result.metadata, step)
             else:
                 self._eval_count += len(values)
 
diff --git a/qiskit_algorithms/gradients/base/base_estimator_gradient.py b/qiskit_algorithms/gradients/base/base_estimator_gradient.py
index f7ea927b..bfc79836 100644
--- a/qiskit_algorithms/gradients/base/base_estimator_gradient.py
+++ b/qiskit_algorithms/gradients/base/base_estimator_gradient.py
@@ -149,7 +149,7 @@ def run(
         job = AlgorithmJob(
             self._run, circuits, observables, parameter_values, parameters, **opts.__dict__
         )
-        job.submit()
+        job._submit()
         return job
 
     @abstractmethod
diff --git a/qiskit_algorithms/gradients/base/base_qgt.py b/qiskit_algorithms/gradients/base/base_qgt.py
index 2e254a8f..24e2e2d1 100644
--- a/qiskit_algorithms/gradients/base/base_qgt.py
+++ b/qiskit_algorithms/gradients/base/base_qgt.py
@@ -161,7 +161,7 @@ def run(
         opts = copy(self._default_options)
         opts.update_options(**options)
         job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **opts.__dict__)
-        job.submit()
+        job._submit()
         return job
 
     @abstractmethod
diff --git a/qiskit_algorithms/gradients/base/base_sampler_gradient.py b/qiskit_algorithms/gradients/base/base_sampler_gradient.py
index 1114b5f0..f17da118 100644
--- a/qiskit_algorithms/gradients/base/base_sampler_gradient.py
+++ b/qiskit_algorithms/gradients/base/base_sampler_gradient.py
@@ -107,7 +107,7 @@ def run(
         opts = copy(self._default_options)
         opts.update_options(**options)
         job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **opts.__dict__)
-        job.submit()
+        job._submit()
         return job
 
     @abstractmethod
diff --git a/qiskit_algorithms/gradients/qfi.py b/qiskit_algorithms/gradients/qfi.py
index 4a53b20d..4bb5623a 100644
--- a/qiskit_algorithms/gradients/qfi.py
+++ b/qiskit_algorithms/gradients/qfi.py
@@ -103,7 +103,7 @@ def run(
         opts = copy(self._default_options)
         opts.update_options(**options)
         job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **opts.__dict__)
-        job.submit()
+        job._submit()
         return job
 
     def _run(
diff --git a/qiskit_algorithms/minimum_eigensolvers/adapt_vqe.py b/qiskit_algorithms/minimum_eigensolvers/adapt_vqe.py
index 895a8fae..bdd22936 100644
--- a/qiskit_algorithms/minimum_eigensolvers/adapt_vqe.py
+++ b/qiskit_algorithms/minimum_eigensolvers/adapt_vqe.py
@@ -62,7 +62,7 @@ class AdaptVQE(VariationalAlgorithm, MinimumEigensolver):
 
       from qiskit_algorithms.minimum_eigensolvers import AdaptVQE, VQE
       from qiskit_algorithms.optimizers import SLSQP
-      from qiskit.primitives import Estimator
+      from qiskit.primitives import StatevectorEstimator
       from qiskit.circuit.library import EvolvedOperatorAnsatz
 
       # get your Hamiltonian
@@ -71,7 +71,7 @@ class AdaptVQE(VariationalAlgorithm, MinimumEigensolver):
       # construct your ansatz
       ansatz = EvolvedOperatorAnsatz(...)
 
-      vqe = VQE(Estimator(), ansatz, SLSQP())
+      vqe = VQE(StatevectorEstimator(), ansatz, SLSQP())
 
       adapt_vqe = AdaptVQE(vqe)
 
diff --git a/qiskit_algorithms/minimum_eigensolvers/diagonal_estimator.py b/qiskit_algorithms/minimum_eigensolvers/diagonal_estimator.py
index aaabf3b8..80b7ef62 100644
--- a/qiskit_algorithms/minimum_eigensolvers/diagonal_estimator.py
+++ b/qiskit_algorithms/minimum_eigensolvers/diagonal_estimator.py
@@ -21,7 +21,7 @@
 
 import numpy as np
 from qiskit.circuit import QuantumCircuit
-from qiskit.primitives import BaseSampler, BaseEstimator, EstimatorResult
+from qiskit.primitives import BaseSamplerV2, BaseEstimatorV2, PubResult
 from qiskit.primitives.utils import init_observable, _circuit_key
 from qiskit.quantum_info import SparsePauliOp
 from qiskit.quantum_info.operators.base_operator import BaseOperator
@@ -30,19 +30,19 @@
 
 
 @dataclass(frozen=True)
-class _DiagonalEstimatorResult(EstimatorResult):
+class _DiagonalEstimatorResult(PubResult):
     """A result from an expectation of a diagonal observable."""
 
     # TODO make each measurement a dataclass rather than a dict
     best_measurements: Sequence[Mapping[str, Any]] | None = None
 
 
-class _DiagonalEstimator(BaseEstimator):
+class _DiagonalEstimator(BaseEstimatorV2):
     """An estimator for diagonal observables."""
 
     def __init__(
         self,
-        sampler: BaseSampler,
+        sampler: BaseSamplerV2,
         aggregation: float | Callable[[Iterable[tuple[float, float]]], float] | None = None,
         callback: Callable[[Sequence[Mapping[str, Any]]], None] | None = None,
         **options,
@@ -104,7 +104,7 @@ def _run(
         job = AlgorithmJob(
             self._call, circuit_indices, observable_indices, parameter_values, **run_options
         )
-        job.submit()
+        job._submit()
         return job
 
     def _call(
@@ -114,13 +114,22 @@ def _call(
         parameter_values: Sequence[Sequence[float]],
         **run_options,
     ) -> _DiagonalEstimatorResult:
+
+        pubs = list(zip(circuits, parameter_values))
         job = self.sampler.run(
-            [self._circuits[i] for i in circuits],
-            parameter_values,
+            pubs,
             **run_options,
         )
-        sampler_result = job.result()
-        samples = sampler_result.quasi_dists
+        results = job.result()
+        samples = []
+        for i, pubres in enumerate(results):
+            qc = pubs[i][0]
+            sampler_result = getattr(results[i].data, qc.cregs[0].name)
+            sample = {
+                label: value / sampler_result.num_shots
+                for label, value in sampler_result.get_counts().items()
+            }
+            samples.append(sample)
 
         # a list of dictionaries containing: {state: (measurement probability, value)}
         evaluations: list[dict[int, tuple[float, float]]] = [
diff --git a/qiskit_algorithms/minimum_eigensolvers/qaoa.py b/qiskit_algorithms/minimum_eigensolvers/qaoa.py
index 8953b795..fbb7a4e6 100644
--- a/qiskit_algorithms/minimum_eigensolvers/qaoa.py
+++ b/qiskit_algorithms/minimum_eigensolvers/qaoa.py
@@ -20,12 +20,13 @@
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library.n_local.qaoa_ansatz import QAOAAnsatz
 from qiskit.quantum_info.operators.base_operator import BaseOperator
-from qiskit.primitives import BaseSampler
+from qiskit.primitives import BaseSamplerV2
 
 from qiskit_algorithms.utils.validation import validate_min
 from qiskit_algorithms.optimizers import Minimizer, Optimizer
 
 from .sampling_vqe import SamplingVQE
+from ..custom_types import Transpiler
 
 
 class QAOA(SamplingVQE):
@@ -54,7 +55,7 @@ class QAOA(SamplingVQE):
     the QAOA object has been constructed.
 
     Attributes:
-        sampler (BaseSampler): The sampler primitive to sample the circuits.
+        sampler (BaseSamplerV2): The sampler primitive to sample the circuits.
         optimizer (Optimizer | Minimizer): A classical optimizer to find the minimum energy. This
             can either be an :class:`.Optimizer` or a callable implementing the
             :class:`.Minimizer` protocol.
@@ -71,6 +72,11 @@ class QAOA(SamplingVQE):
             that can access the intermediate data at each optimization step. These data are: the
             evaluation count, the optimizer parameters for the ansatz, the evaluated value, and
             the metadata dictionary.
+        transpiler: An optional object with a `run` method allowing to transpile the circuits
+            that are produced within this algorithm. If set to `None`, these won't be transpiled.
+        transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+            method as keyword arguments.
+
 
     References:
         [1]: Farhi, E., Goldstone, J., Gutmann, S., "A Quantum Approximate Optimization Algorithm"
@@ -85,7 +91,7 @@ class QAOA(SamplingVQE):
 
     def __init__(
         self,
-        sampler: BaseSampler,
+        sampler: BaseSamplerV2,
         optimizer: Optimizer | Minimizer,
         *,
         reps: int = 1,
@@ -94,6 +100,8 @@ def __init__(
         initial_point: np.ndarray | None = None,
         aggregation: float | Callable[[list[float]], float] | None = None,
         callback: Callable[[int, np.ndarray, float, dict[str, Any]], None] | None = None,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
     ) -> None:
         r"""
         Args:
@@ -117,6 +125,10 @@ def __init__(
             callback: A callback that can access the intermediate data at each optimization step.
                 These data are: the evaluation count, the optimizer parameters for the ansatz, the
                 evaluated value, the metadata dictionary.
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced within this algorithm. If set to `None`, these won't be transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
         """
         validate_min("reps", reps, 1)
 
@@ -124,6 +136,8 @@ def __init__(
         self.mixer = mixer
         self.initial_state = initial_state
         self._cost_operator = None
+        self._transpiler = transpiler
+        self._transpiler_options = transpiler_options if transpiler_options is not None else {}
 
         super().__init__(
             sampler=sampler,
@@ -136,6 +150,9 @@ def __init__(
 
     def _check_operator_ansatz(self, operator: BaseOperator):
         # Recreates a circuit based on operator parameter.
-        self.ansatz = QAOAAnsatz(
+        ansatz = QAOAAnsatz(
             operator, self.reps, initial_state=self.initial_state, mixer_operator=self.mixer
         ).decompose()  # TODO remove decompose once #6674 is fixed
+        if self._transpiler is not None:
+            ansatz = self._transpiler.run(ansatz, **self._transpiler_options)
+        self.ansatz = ansatz
diff --git a/qiskit_algorithms/minimum_eigensolvers/sampling_vqe.py b/qiskit_algorithms/minimum_eigensolvers/sampling_vqe.py
index 0ebd7e05..ae25976b 100644
--- a/qiskit_algorithms/minimum_eigensolvers/sampling_vqe.py
+++ b/qiskit_algorithms/minimum_eigensolvers/sampling_vqe.py
@@ -22,7 +22,7 @@
 import numpy as np
 
 from qiskit.circuit import QuantumCircuit
-from qiskit.primitives import BaseSampler
+from qiskit.primitives import BaseSamplerV2
 from qiskit.result import QuasiDistribution
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 
@@ -92,7 +92,7 @@ def my_minimizer(fun, x0, jac=None, bounds=None) -> OptimizerResult:
     the ``SamplingVQE`` object has been constructed.
 
     Attributes:
-        sampler (BaseSampler): The sampler primitive to sample the circuits.
+        sampler (BaseSamplerV2): The sampler primitive to sample the circuits.
         ansatz (QuantumCircuit): A parameterized quantum circuit to prepare the trial state.
         optimizer (Optimizer | Minimizer): A classical optimizer to find the minimum energy. This
             can either be an :class:`.Optimizer` or a callable implementing the
@@ -116,7 +116,7 @@ def my_minimizer(fun, x0, jac=None, bounds=None) -> OptimizerResult:
 
     def __init__(
         self,
-        sampler: BaseSampler,
+        sampler: BaseSamplerV2,
         ansatz: QuantumCircuit,
         optimizer: Optimizer | Minimizer,
         *,
@@ -240,7 +240,12 @@ def compute_minimum_eigenvalue(
             optimizer_result.x,
         )
 
-        final_state = self.sampler.run([self.ansatz], [optimizer_result.x]).result().quasi_dists[0]
+        final_res = self.sampler.run([(self.ansatz, optimizer_result.x)]).result()
+        final_state = getattr(final_res[0].data, self.ansatz.cregs[0].name)
+        final_state = {
+            label: value / final_state.num_shots
+            for label, value in final_state.get_counts().items()
+        }
 
         if aux_operators is not None:
             aux_operators_evaluated = estimate_observables(
@@ -313,13 +318,14 @@ def store_best_measurement(best):
         def evaluate_energy(parameters: np.ndarray) -> np.ndarray | float:
             nonlocal eval_count
             # handle broadcasting: ensure parameters is of shape [array, array, ...]
-            parameters = np.reshape(parameters, (-1, num_parameters)).tolist()
-            batch_size = len(parameters)
-
-            estimator_result = estimator.run(
-                batch_size * [ansatz], batch_size * [operator], parameters
-            ).result()
-            values = estimator_result.values
+            # parameters = np.reshape(parameters, (-1, num_parameters)).tolist()
+            # batch_size = len(parameters)
+
+            job = self.estimator.run([(ansatz, operator, parameters)])
+            estimator_result = job.result()[0]
+            values = estimator_result.data.evs
+            if not values.shape:
+                values = values.reshape(1)
 
             if self.callback is not None:
                 metadata = estimator_result.metadata
diff --git a/qiskit_algorithms/minimum_eigensolvers/vqe.py b/qiskit_algorithms/minimum_eigensolvers/vqe.py
index b0e85a67..c449da4f 100644
--- a/qiskit_algorithms/minimum_eigensolvers/vqe.py
+++ b/qiskit_algorithms/minimum_eigensolvers/vqe.py
@@ -22,7 +22,7 @@
 import numpy as np
 
 from qiskit.circuit import QuantumCircuit
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 
 from qiskit_algorithms.gradients import BaseEstimatorGradient
@@ -94,7 +94,7 @@ def my_minimizer(fun, x0, jac=None, bounds=None) -> OptimizerResult:
     the VQE object has been constructed.
 
     Attributes:
-        estimator (BaseEstimator): The estimator primitive to compute the expectation value of the
+        estimator (BaseEstimatorV2): The estimator primitive to compute the expectation value of the
             Hamiltonian operator.
         ansatz (QuantumCircuit): A parameterized quantum circuit to prepare the trial state.
         optimizer (Optimizer | Minimizer): A classical optimizer to find the minimum energy. This
@@ -114,7 +114,7 @@ def my_minimizer(fun, x0, jac=None, bounds=None) -> OptimizerResult:
 
     def __init__(
         self,
-        estimator: BaseEstimator,
+        estimator: BaseEstimatorV2,
         ansatz: QuantumCircuit,
         optimizer: Optimizer | Minimizer,
         *,
@@ -258,16 +258,19 @@ def evaluate_energy(parameters: np.ndarray) -> np.ndarray | float:
             nonlocal eval_count
 
             # handle broadcasting: ensure parameters is of shape [array, array, ...]
-            parameters = np.reshape(parameters, (-1, num_parameters)).tolist()
-            batch_size = len(parameters)
+            # parameters = np.reshape(parameters, (-1, num_parameters)).tolist()
+            # batch_size = len(parameters)
 
             try:
-                job = self.estimator.run(batch_size * [ansatz], batch_size * [operator], parameters)
-                estimator_result = job.result()
+                job = self.estimator.run([(ansatz, operator, parameters)])
+                estimator_result = job.result()[0]
             except Exception as exc:
                 raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc
 
-            values = estimator_result.values
+            values = estimator_result.data.evs
+
+            if not values.shape:
+                values = values.reshape(1)
 
             if self.callback is not None:
                 metadata = estimator_result.metadata
diff --git a/qiskit_algorithms/observables_evaluator.py b/qiskit_algorithms/observables_evaluator.py
index ae125bfb..a455c78e 100644
--- a/qiskit_algorithms/observables_evaluator.py
+++ b/qiskit_algorithms/observables_evaluator.py
@@ -20,7 +20,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.quantum_info import SparsePauliOp
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 
 from .exceptions import AlgorithmError
@@ -28,7 +28,7 @@
 
 
 def estimate_observables(
-    estimator: BaseEstimator,
+    estimator: BaseEstimatorV2,
     quantum_state: QuantumCircuit,
     observables: ListOrDict[BaseOperator],
     parameter_values: Sequence[float] | None = None,
@@ -42,7 +42,7 @@ def estimate_observables(
     Args:
         estimator: An estimator primitive used for calculations.
         quantum_state: A (parameterized) quantum circuit preparing a quantum state that expectation
-            values are computed against.
+            values are computed against. It is expected to be an ISA circuit.
         observables: A list or a dictionary of operators whose expectation values are to be
             calculated.
         parameter_values: Optional list of parameters values to evaluate the quantum circuit on.
@@ -63,21 +63,19 @@ def estimate_observables(
 
     if len(observables_list) > 0:
         observables_list = _handle_zero_ops(observables_list)
-        quantum_state = [quantum_state] * len(observables)
         parameter_values_: Sequence[float] | Sequence[Sequence[float]] | None = parameter_values
-        if parameter_values is not None:
-            parameter_values_ = [parameter_values] * len(observables)
         try:
-            estimator_job = estimator.run(quantum_state, observables_list, parameter_values_)
-            expectation_values = estimator_job.result().values
+            estimator_job = estimator.run([(quantum_state, observables_list, parameter_values_)])
+            estimator_result = estimator_job.result()[0]
+            expectation_values = estimator_result.data.evs
         except Exception as exc:
             raise AlgorithmError("The primitive job failed!") from exc
 
-        metadata = estimator_job.result().metadata
+        metadata = estimator_result.metadata
         # Discard values below threshold
         observables_means = expectation_values * (np.abs(expectation_values) > threshold)
         # zip means and metadata into tuples
-        observables_results = list(zip(observables_means, metadata))
+        observables_results = list(zip(observables_means, [metadata] * len(observables_means)))
     else:
         observables_results = []
 
diff --git a/qiskit_algorithms/optimizers/qnspsa.py b/qiskit_algorithms/optimizers/qnspsa.py
index 4c00741d..45659594 100644
--- a/qiskit_algorithms/optimizers/qnspsa.py
+++ b/qiskit_algorithms/optimizers/qnspsa.py
@@ -20,10 +20,11 @@
 import numpy as np
 from qiskit.circuit import QuantumCircuit
 
-from qiskit.primitives import BaseSampler
+from qiskit.primitives import BaseSamplerV2
 from qiskit_algorithms.state_fidelities import ComputeUncompute
 
 from .spsa import SPSA, CALLBACK, TERMINATIONCHECKER, _batch_evaluate
+from ..custom_types import Transpiler
 
 # the function to compute the fidelity
 FIDELITY = Callable[[np.ndarray, np.ndarray], float]
@@ -63,7 +64,8 @@ class QNSPSA(SPSA):
             import numpy as np
             from qiskit_algorithms.optimizers import QNSPSA
             from qiskit.circuit.library import PauliTwoDesign
-            from qiskit.primitives import Estimator, Sampler
+            from qiskit.primitives import StatevectorEstimator as Estimator,\
+                StatevectorSampler as Sampler
             from qiskit.quantum_info import Pauli
 
             # problem setup
@@ -75,8 +77,8 @@ class QNSPSA(SPSA):
             estimator = Estimator()
 
             def loss(x):
-                result = estimator.run([ansatz], [observable], [x]).result()
-                return np.real(result.values[0])
+                result = estimator.run([(ansatz, observable, x)]).result()[0]
+                return np.real(result.data.evs[0])
 
             # fidelity for estimation of the geometric tensor
             sampler = Sampler()
@@ -233,7 +235,9 @@ def settings(self) -> dict[str, Any]:
     def get_fidelity(
         circuit: QuantumCircuit,
         *,
-        sampler: BaseSampler | None = None,
+        sampler: BaseSamplerV2 | None = None,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
     ) -> Callable[[np.ndarray, np.ndarray], float]:
         r"""Get a function to compute the fidelity of ``circuit`` with itself.
 
@@ -251,12 +255,19 @@ def get_fidelity(
         Args:
             circuit: The circuit preparing the parameterized ansatz.
             sampler: A sampler primitive to sample from a quantum state.
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced by the fidelity object. If set to `None`, these won't be
+                transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
 
         Returns:
             A handle to the function :math:`F`.
 
         """
-        fid = ComputeUncompute(sampler)
+        fid = ComputeUncompute(
+            sampler, transpiler=transpiler, transpiler_options=transpiler_options
+        )
 
         num_parameters = circuit.num_parameters
 
diff --git a/qiskit_algorithms/optimizers/spsa.py b/qiskit_algorithms/optimizers/spsa.py
index 1a13fb97..1b0633ff 100644
--- a/qiskit_algorithms/optimizers/spsa.py
+++ b/qiskit_algorithms/optimizers/spsa.py
@@ -90,7 +90,7 @@ class SPSA(Optimizer):
             import numpy as np
             from qiskit_algorithms.optimizers import SPSA
             from qiskit.circuit.library import PauliTwoDesign
-            from qiskit.primitives import Estimator
+            from qiskit.primitives import StatevectorEstimator as Estimator
             from qiskit.quantum_info import SparsePauliOp
 
             ansatz = PauliTwoDesign(2, reps=1, seed=2)
@@ -99,8 +99,8 @@ class SPSA(Optimizer):
             estimator = Estimator()
 
             def loss(x):
-                job = estimator.run([ansatz], [observable], [x])
-                return job.result().values[0]
+                job = estimator.run([(ansatz, observable, x)])
+                return job.result()[0].data.evs
 
             spsa = SPSA(maxiter=300)
             result = spsa.minimize(loss, x0=initial_point)
diff --git a/qiskit_algorithms/optimizers/umda.py b/qiskit_algorithms/optimizers/umda.py
index f420e82d..2ad7967d 100644
--- a/qiskit_algorithms/optimizers/umda.py
+++ b/qiskit_algorithms/optimizers/umda.py
@@ -74,7 +74,7 @@ class UMDA(Optimizer):
             from qiskit_algorithms.optimizers import UMDA
             from qiskit_algorithms import QAOA
             from qiskit.quantum_info import Pauli
-            from qiskit.primitives import Sampler
+            from qiskit.primitives import StatevectorSampler as Sampler
 
             X = Pauli("X")
             I = Pauli("I")
diff --git a/qiskit_algorithms/phase_estimators/hamiltonian_phase_estimation.py b/qiskit_algorithms/phase_estimators/hamiltonian_phase_estimation.py
index bfb36e9a..35711762 100644
--- a/qiskit_algorithms/phase_estimators/hamiltonian_phase_estimation.py
+++ b/qiskit_algorithms/phase_estimators/hamiltonian_phase_estimation.py
@@ -1,6 +1,6 @@
 # This code is part of a Qiskit project.
 #
-# (C) Copyright IBM 2020, 2023.
+# (C) Copyright IBM 2020, 2024.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -14,16 +14,18 @@
 
 from __future__ import annotations
 
+from typing import Any
 
 from qiskit import QuantumCircuit
 from qiskit.circuit.library import PauliEvolutionGate
-from qiskit.primitives import BaseSampler
+from qiskit.primitives import BaseSamplerV2
 from qiskit.quantum_info import SparsePauliOp, Statevector, Pauli
 from qiskit.synthesis import EvolutionSynthesis
 
-from .phase_estimation import PhaseEstimation
 from .hamiltonian_phase_estimation_result import HamiltonianPhaseEstimationResult
+from .phase_estimation import PhaseEstimation
 from .phase_estimation_scale import PhaseEstimationScale
+from ..custom_types import Transpiler
 
 
 class HamiltonianPhaseEstimation:
@@ -83,17 +85,26 @@ class HamiltonianPhaseEstimation:
     def __init__(
         self,
         num_evaluation_qubits: int,
-        sampler: BaseSampler | None = None,
+        sampler: BaseSamplerV2 | None = None,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
     ) -> None:
         r"""
         Args:
             num_evaluation_qubits: The number of qubits used in estimating the phase. The phase will
                 be estimated as a binary string with this many bits.
             sampler: The sampler primitive on which the circuit will be sampled.
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced within this algorithm. If set to `None`, these won't be
+                transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
         """
         self._phase_estimation = PhaseEstimation(
             num_evaluation_qubits=num_evaluation_qubits,
             sampler=sampler,
+            transpiler=transpiler,
+            transpiler_options=transpiler_options,
         )
 
     def _get_scale(self, hamiltonian, bound=None) -> PhaseEstimationScale:
diff --git a/qiskit_algorithms/phase_estimators/ipe.py b/qiskit_algorithms/phase_estimators/ipe.py
index bac41756..991909e0 100644
--- a/qiskit_algorithms/phase_estimators/ipe.py
+++ b/qiskit_algorithms/phase_estimators/ipe.py
@@ -15,16 +15,19 @@
 
 from __future__ import annotations
 
+from typing import Any
+
 import numpy
 
 from qiskit.circuit import QuantumCircuit, QuantumRegister
 from qiskit.circuit.classicalregister import ClassicalRegister
-from qiskit.primitives import BaseSampler
+from qiskit.primitives import BaseSamplerV2
 
 from qiskit_algorithms.exceptions import AlgorithmError
 
 from .phase_estimator import PhaseEstimator
 from .phase_estimator import PhaseEstimatorResult
+from ..custom_types import Transpiler
 
 
 class IterativePhaseEstimation(PhaseEstimator):
@@ -40,12 +43,19 @@ class IterativePhaseEstimation(PhaseEstimator):
     def __init__(
         self,
         num_iterations: int,
-        sampler: BaseSampler | None = None,
+        sampler: BaseSamplerV2 | None = None,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
     ) -> None:
         r"""
         Args:
             num_iterations: The number of iterations (rounds) of the phase estimation to run.
             sampler: The sampler primitive on which the circuit will be sampled.
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced within this algorithm. If set to `None`, these won't be
+                transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
 
         Raises:
             ValueError: if num_iterations is not greater than zero.
@@ -58,6 +68,8 @@ def __init__(
             raise ValueError("`num_iterations` must be greater than zero.")
         self._num_iterations = num_iterations
         self._sampler = sampler
+        self._transpiler = transpiler
+        self._transpiler_options = transpiler_options if transpiler_options is not None else {}
 
     # pylint: disable=too-many-positional-arguments
     def construct_circuit(
@@ -126,9 +138,17 @@ def _estimate_phase_iteratively(self, unitary, state_preparation):
             qc = self.construct_circuit(
                 unitary, state_preparation, k, -2 * numpy.pi * omega_coef, True
             )
+
+            if self._transpiler is not None:
+                qc = self._transpiler.run(qc, **self._transpiler_options)
+
             try:
                 sampler_job = self._sampler.run([qc])
-                result = sampler_job.result().quasi_dists[0]
+                result = sampler_job.result()[0].data.c
+                result = {
+                    label: value / result.num_shots
+                    for label, value in result.get_int_counts().items()
+                }
             except Exception as exc:
                 raise AlgorithmError("The primitive job failed!") from exc
             x = 1 if result.get(1, 0) > result.get(0, 0) else 0
diff --git a/qiskit_algorithms/phase_estimators/phase_estimation.py b/qiskit_algorithms/phase_estimators/phase_estimation.py
index bf84b736..32cc2812 100644
--- a/qiskit_algorithms/phase_estimators/phase_estimation.py
+++ b/qiskit_algorithms/phase_estimators/phase_estimation.py
@@ -1,6 +1,6 @@
 # This code is part of a Qiskit project.
 #
-# (C) Copyright IBM 2020, 2023.
+# (C) Copyright IBM 2020, 2024.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -15,18 +15,21 @@
 
 from __future__ import annotations
 
+from typing import Any
+
 import numpy
 import qiskit
 from qiskit import circuit
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.classicalregister import ClassicalRegister
-from qiskit.primitives import BaseSampler
+from qiskit.primitives import BaseSamplerV2
 from qiskit.result import Result
 
 from qiskit_algorithms.exceptions import AlgorithmError
 
 from .phase_estimation_result import PhaseEstimationResult, _sort_phases
 from .phase_estimator import PhaseEstimator
+from ..custom_types import Transpiler
 
 
 class PhaseEstimation(PhaseEstimator):
@@ -82,13 +85,20 @@ class PhaseEstimation(PhaseEstimator):
     def __init__(
         self,
         num_evaluation_qubits: int,
-        sampler: BaseSampler | None = None,
+        sampler: BaseSamplerV2 | None = None,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
     ) -> None:
         r"""
         Args:
             num_evaluation_qubits: The number of qubits used in estimating the phase. The phase will
                 be estimated as a binary string with this many bits.
             sampler: The sampler primitive on which the circuit will be sampled.
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced within this algorithm. If set to `None`, these won't be
+                transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
 
         Raises:
             AlgorithmError: If a sampler is not provided
@@ -101,6 +111,8 @@ def __init__(
             self._num_evaluation_qubits = num_evaluation_qubits
 
         self._sampler = sampler
+        self._transpiler = transpiler
+        self._transpiler_options = transpiler_options if transpiler_options is not None else {}
 
     def construct_circuit(
         self, unitary: QuantumCircuit, state_preparation: QuantumCircuit | None = None
@@ -189,6 +201,9 @@ def estimate_from_pe_circuit(self, pe_circuit: QuantumCircuit) -> PhaseEstimatio
             AlgorithmError: Primitive job failed.
         """
 
+        if self._transpiler is not None:
+            pe_circuit = self._transpiler.run(pe_circuit, **self._transpiler_options)
+
         self._add_measurement_if_required(pe_circuit)
 
         try:
@@ -196,11 +211,13 @@ def estimate_from_pe_circuit(self, pe_circuit: QuantumCircuit) -> PhaseEstimatio
             circuit_result = circuit_job.result()
         except Exception as exc:
             raise AlgorithmError("The primitive job failed!") from exc
-        phases = circuit_result.quasi_dists[0]
+        phases = circuit_result[0].data.meas.get_counts()
+        # Ensure we still return the measurement strings in sorted order, which SamplerV2 doesn't
+        # guarantee
+        measurement_labels = sorted(phases.keys())
         phases_bitstrings = {}
-        for key, phase in phases.items():
-            bitstring_key = self._get_reversed_bitstring(self._num_evaluation_qubits, key)
-            phases_bitstrings[bitstring_key] = phase
+        for key in measurement_labels:
+            phases_bitstrings[key[::-1]] = phases[key] / circuit_result[0].data.meas.num_shots
         phases = phases_bitstrings
 
         return PhaseEstimationResult(
diff --git a/qiskit_algorithms/phase_estimators/phase_estimation_result.py b/qiskit_algorithms/phase_estimators/phase_estimation_result.py
index bf57c800..31df327c 100644
--- a/qiskit_algorithms/phase_estimators/phase_estimation_result.py
+++ b/qiskit_algorithms/phase_estimators/phase_estimation_result.py
@@ -23,7 +23,7 @@ class PhaseEstimationResult(PhaseEstimatorResult):
 
     This class is instantiated by the ``PhaseEstimation`` class, not via user code.
     The ``PhaseEstimation`` class generates a list of phases and corresponding weights. Upon
-    completion it returns the results as an instance of this class. The main method for
+    completion, it returns the results as an instance of this class. The main method for
     accessing the results is `filter_phases`.
 
     The canonical phase satisfying the ``PhaseEstimator`` interface, returned by the
diff --git a/qiskit_algorithms/phase_estimators/phase_estimator.py b/qiskit_algorithms/phase_estimators/phase_estimator.py
index 1f0f5002..2a78f7f8 100644
--- a/qiskit_algorithms/phase_estimators/phase_estimator.py
+++ b/qiskit_algorithms/phase_estimators/phase_estimator.py
@@ -1,6 +1,6 @@
 # This code is part of a Qiskit project.
 #
-# (C) Copyright IBM 2020, 2023.
+# (C) Copyright IBM 2020, 2024.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -38,10 +38,6 @@ def estimate(
         """Estimate the phase."""
         raise NotImplementedError
 
-    @staticmethod
-    def _get_reversed_bitstring(length: int, number: int) -> str:
-        return f"{number:b}".zfill(length)[::-1]
-
 
 class PhaseEstimatorResult(AlgorithmResult):
     """Phase Estimator Result."""
diff --git a/qiskit_algorithms/state_fidelities/base_state_fidelity.py b/qiskit_algorithms/state_fidelities/base_state_fidelity.py
index 5c1199c4..d6a0a864 100644
--- a/qiskit_algorithms/state_fidelities/base_state_fidelity.py
+++ b/qiskit_algorithms/state_fidelities/base_state_fidelity.py
@@ -16,7 +16,7 @@
 from __future__ import annotations
 from abc import ABC, abstractmethod
 from collections.abc import MutableMapping
-from typing import cast, Sequence, List
+from typing import cast, Sequence, List, Any
 import numpy as np
 
 from qiskit import QuantumCircuit
@@ -24,6 +24,7 @@
 from qiskit.primitives.utils import _circuit_key
 
 from ..algorithm_job import AlgorithmJob
+from ..custom_types import Transpiler
 
 
 class BaseStateFidelity(ABC):
@@ -42,10 +43,23 @@ class BaseStateFidelity(ABC):
 
     """
 
-    def __init__(self) -> None:
-
+    def __init__(
+        self,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
+    ) -> None:
+        r"""
+        Args:
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced within this algorithm. If set to `None`, these won't be
+                transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
+        """
         # use cache for preventing unnecessary circuit compositions
         self._circuit_cache: MutableMapping[tuple[int, int], QuantumCircuit] = {}
+        self._transpiler = transpiler
+        self._transpiler_options = transpiler_options if transpiler_options is not None else {}
 
     @staticmethod
     def _preprocess_values(
@@ -195,6 +209,9 @@ def _construct_circuits(
                 # update cache
                 self._circuit_cache[_circuit_key(circuit_1), _circuit_key(circuit_2)] = circuit
 
+        if self._transpiler is not None:
+            return self._transpiler.run(circuits, **self._transpiler_options)
+
         return circuits
 
     def _construct_value_list(
@@ -245,7 +262,7 @@ def _run(
         circuits_2: QuantumCircuit | Sequence[QuantumCircuit],
         values_1: Sequence[float] | Sequence[Sequence[float]] | None = None,
         values_2: Sequence[float] | Sequence[Sequence[float]] | None = None,
-        **options,
+        shots: int | Sequence[int] | None = None,
     ) -> AlgorithmJob:
         r"""
         Computes the state overlap (fidelity) calculation between two
@@ -257,10 +274,12 @@ def _run(
             circuits_2: (Parametrized) quantum circuits preparing :math:`|\phi\rangle`.
             values_1: Numerical parameters to be bound to the first set of circuits
             values_2: Numerical parameters to be bound to the second set of circuits.
-            options: Primitive backend runtime options used for circuit execution. The order
-                of priority is\: options in ``run`` method > fidelity's default
-                options > primitive's default setting.
-                Higher priority setting overrides lower priority setting.
+            shots: Number of shots to be used by the underlying sampler. If a single integer is
+                provided, this number will be used for all circuits. If a sequence of integers is
+                provided, they will be used on a per-circuit basis. If none is provided, the
+                fidelity's default number of shots will be used for all circuits. If this number is
+                also set to None, the underlying primitive's default number of shots will be used
+                for all circuits.
 
         Returns:
             A newly constructed algorithm job instance to get the fidelity result.
@@ -273,7 +292,7 @@ def run(
         circuits_2: QuantumCircuit | Sequence[QuantumCircuit],
         values_1: Sequence[float] | Sequence[Sequence[float]] | None = None,
         values_2: Sequence[float] | Sequence[Sequence[float]] | None = None,
-        **options,
+        shots: int | Sequence[int] | None = None,
     ) -> AlgorithmJob:
         r"""
         Runs asynchronously the state overlap (fidelity) calculation between two
@@ -286,18 +305,20 @@ def run(
             circuits_2: (Parametrized) quantum circuits preparing :math:`|\phi\rangle`.
             values_1: Numerical parameters to be bound to the first set of circuits.
             values_2: Numerical parameters to be bound to the second set of circuits.
-            options: Primitive backend runtime options used for circuit execution. The order
-                of priority is\: options in ``run`` method > fidelity's default
-                options > primitive's default setting.
-                Higher priority setting overrides lower priority setting.
+            shots: Number of shots to be used by the underlying sampler. If a single integer is
+                provided, this number will be used for all circuits. If a sequence of integers is
+                provided, they will be used on a per-circuit basis. If none is provided, the
+                fidelity's default number of shots will be used for all circuits. If this number is
+                also set to None, the underlying primitive's default number of shots will be used
+                for all circuits.
 
         Returns:
             Primitive job for the fidelity calculation.
             The job's result is an instance of :class:`.StateFidelityResult`.
         """
-        job = self._run(circuits_1, circuits_2, values_1, values_2, **options)
+        job = self._run(circuits_1, circuits_2, values_1, values_2, shots)
 
-        job.submit()
+        job._submit()
         return job
 
     @staticmethod
diff --git a/qiskit_algorithms/state_fidelities/compute_uncompute.py b/qiskit_algorithms/state_fidelities/compute_uncompute.py
index b16e4011..8bdba579 100644
--- a/qiskit_algorithms/state_fidelities/compute_uncompute.py
+++ b/qiskit_algorithms/state_fidelities/compute_uncompute.py
@@ -14,18 +14,20 @@
 """
 
 from __future__ import annotations
+
 from collections.abc import Sequence
-from copy import copy
+from typing import Any
 
 from qiskit import QuantumCircuit
-from qiskit.primitives import BaseSampler
+from qiskit.primitives import BaseSamplerV2
+from qiskit.primitives.containers.sampler_pub import SamplerPub
 from qiskit.primitives.primitive_job import PrimitiveJob
-from qiskit.providers import Options
 
-from ..exceptions import AlgorithmError
 from .base_state_fidelity import BaseStateFidelity
 from .state_fidelity_result import StateFidelityResult
 from ..algorithm_job import AlgorithmJob
+from ..custom_types import Transpiler
+from ..exceptions import AlgorithmError
 
 
 class ComputeUncompute(BaseStateFidelity):
@@ -53,16 +55,18 @@ class ComputeUncompute(BaseStateFidelity):
 
     def __init__(
         self,
-        sampler: BaseSampler,
-        options: Options | None = None,
+        sampler: BaseSamplerV2,
+        shots: int | None = None,
         local: bool = False,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
     ) -> None:
         r"""
         Args:
             sampler: Sampler primitive instance.
-            options: Primitive backend runtime options used for circuit execution.
-                The order of priority is: options in ``run`` method > fidelity's
-                default options > primitive's default setting.
+            shots: Number of shots to be used by the underlying sampler.
+                The order of priority is: number of shots in ``run`` method > fidelity's
+                number of shots > primitive's default number of shots.
                 Higher priority setting overrides lower priority setting.
             local: If set to ``True``, the fidelity is averaged over
                 single-qubit projectors
@@ -75,20 +79,23 @@ def __init__(
                 This coincides with the standard (global) fidelity in the limit of
                 the fidelity approaching 1. Might be used to increase the variance
                 to improve trainability in algorithms such as :class:`~.time_evolvers.PVQD`.
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced within this algorithm. If set to `None`, these won't be
+                transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
 
         Raises:
-            ValueError: If the sampler is not an instance of ``BaseSampler``.
+            ValueError: If the sampler is not an instance of ``BaseSamplerV2``.
         """
-        if not isinstance(sampler, BaseSampler):
+        if not isinstance(sampler, BaseSamplerV2):
             raise ValueError(
-                f"The sampler should be an instance of BaseSampler, " f"but got {type(sampler)}"
+                f"The sampler should be an instance of BaseSamplerV2, " f"but got {type(sampler)}"
             )
-        self._sampler: BaseSampler = sampler
+        self._sampler: BaseSamplerV2 = sampler
         self._local = local
-        self._default_options = Options()
-        if options is not None:
-            self._default_options.update_options(**options)
-        super().__init__()
+        self._shots = shots
+        super().__init__(transpiler, transpiler_options)
 
     def create_fidelity_circuit(
         self, circuit_1: QuantumCircuit, circuit_2: QuantumCircuit
@@ -119,7 +126,7 @@ def _run(
         circuits_2: QuantumCircuit | Sequence[QuantumCircuit],
         values_1: Sequence[float] | Sequence[Sequence[float]] | None = None,
         values_2: Sequence[float] | Sequence[Sequence[float]] | None = None,
-        **options,
+        shots: int | Sequence[int] | None = None,
     ) -> AlgorithmJob:
         r"""
         Computes the state overlap (fidelity) calculation between two
@@ -131,10 +138,12 @@ def _run(
             circuits_2: (Parametrized) quantum circuits preparing :math:`|\phi\rangle`.
             values_1: Numerical parameters to be bound to the first circuits.
             values_2: Numerical parameters to be bound to the second circuits.
-            options: Primitive backend runtime options used for circuit execution.
-                    The order of priority is: options in ``run`` method > fidelity's
-                    default options > primitive's default setting.
-                    Higher priority setting overrides lower priority setting.
+            shots: Number of shots to be used by the underlying sampler. If a single integer is
+                provided, this number will be used for all circuits. If a sequence of integers is
+                provided, they will be used on a per-circuit basis. If none is provided, the
+                fidelity's default number of shots will be used for all circuits. If this number is
+                also set to None, the underlying primitive's default number of shots will be used
+                for all circuits.
 
         Returns:
             An AlgorithmJob for the fidelity calculation.
@@ -143,7 +152,6 @@ def _run(
             ValueError: At least one pair of circuits must be defined.
             AlgorithmError: If the sampler job is not completed successfully.
         """
-
         circuits = self._construct_circuits(circuits_1, circuits_2)
         if len(circuits) == 0:
             raise ValueError(
@@ -151,79 +159,88 @@ def _run(
             )
         values = self._construct_value_list(circuits_1, circuits_2, values_1, values_2)
 
-        # The priority of run options is as follows:
-        # options in `evaluate` method > fidelity's default options >
-        # primitive's default options.
-        opts = copy(self._default_options)
-        opts.update_options(**options)
+        # The priority of number of shots options is as follows:
+        # number in `run` method > fidelity's default number of shots >
+        # primitive's default number of shots.
+        if not isinstance(shots, Sequence):
+            if shots is None:
+                shots = self.shots
+            coerced_pubs = [
+                SamplerPub.coerce((circuit, value), shots)
+                for circuit, value in zip(circuits, values)
+            ]
+        else:
+            coerced_pubs = [
+                SamplerPub.coerce((circuit, value), shots_number)
+                for circuit, value, shots_number in zip(circuits, values, shots)
+            ]
 
-        sampler_job = self._sampler.run(circuits=circuits, parameter_values=values, **opts.__dict__)
+        job = self._sampler.run(coerced_pubs)
 
-        local_opts = self._get_local_options(opts.__dict__)
-        return AlgorithmJob(ComputeUncompute._call, sampler_job, circuits, self._local, local_opts)
+        return AlgorithmJob(ComputeUncompute._call, job, circuits, self._local)
 
     @staticmethod
     def _call(
-        job: PrimitiveJob, circuits: Sequence[QuantumCircuit], local: bool, local_opts: Options
+        job: PrimitiveJob, circuits: Sequence[QuantumCircuit], local: bool
     ) -> StateFidelityResult:
         try:
             result = job.result()
         except Exception as exc:
             raise AlgorithmError("Sampler job failed!") from exc
 
+        pub_results_data = [
+            getattr(pub_result.data, circuit.cregs[0].name)
+            for pub_result, circuit in zip(result, circuits)
+        ]
+        quasi_dists = [
+            {
+                label: value / prob_dist.num_shots
+                for label, value in prob_dist.get_int_counts().items()
+            }
+            for prob_dist in pub_results_data
+        ]
+
         if local:
             raw_fidelities = [
                 ComputeUncompute._get_local_fidelity(prob_dist, circuit.num_qubits)
-                for prob_dist, circuit in zip(result.quasi_dists, circuits)
+                for prob_dist, circuit in zip(quasi_dists, circuits)
             ]
         else:
             raw_fidelities = [
-                ComputeUncompute._get_global_fidelity(prob_dist) for prob_dist in result.quasi_dists
+                ComputeUncompute._get_global_fidelity(prob_dist) for prob_dist in quasi_dists
             ]
         fidelities = ComputeUncompute._truncate_fidelities(raw_fidelities)
+        shots = [pub_result_data.num_shots for pub_result_data in pub_results_data]
+
+        if len(shots) == 1:
+            shots = shots[0]
 
         return StateFidelityResult(
             fidelities=fidelities,
             raw_fidelities=raw_fidelities,
             metadata=result.metadata,
-            options=local_opts,
+            shots=shots,
         )
 
     @property
-    def options(self) -> Options:
-        """Return the union of estimator options setting and fidelity default options,
-        where, if the same field is set in both, the fidelity's default options override
-        the primitive's default setting.
+    def shots(self) -> int | None:
+        """Return the number of shots used by the `run` method of the Sampler primitive. If None,
+        the default number of shots of the primitive is used.
 
         Returns:
-            The fidelity default + estimator options.
+            The default number of shots.
         """
-        return self._get_local_options(self._default_options.__dict__)
+        return self._shots
 
-    def update_default_options(self, **options):
-        """Update the fidelity's default options setting.
+    @shots.setter
+    def shots(self, shots: int | None):
+        """Update the fidelity's default number of shots setting.
 
         Args:
-            **options: The fields to update the default options.
+            shots: The new default number of shots.
         """
 
-        self._default_options.update_options(**options)
-
-    def _get_local_options(self, options: Options) -> Options:
-        """Return the union of the primitive's default setting,
-        the fidelity default options, and the options in the ``run`` method.
-        The order of priority is: options in ``run`` method > fidelity's
-                default options > primitive's default setting.
-
-        Args:
-            options: The fields to update the options
-
-        Returns:
-            The fidelity default + estimator + run options.
-        """
-        opts = copy(self._sampler.options)
-        opts.update_options(**options)
-        return opts
+        self._shots = shots
 
     @staticmethod
     def _get_global_fidelity(probability_distribution: dict[int, float]) -> float:
diff --git a/qiskit_algorithms/state_fidelities/state_fidelity_result.py b/qiskit_algorithms/state_fidelities/state_fidelity_result.py
index 6dc26dbf..c0f7bd9a 100644
--- a/qiskit_algorithms/state_fidelities/state_fidelity_result.py
+++ b/qiskit_algorithms/state_fidelities/state_fidelity_result.py
@@ -16,10 +16,8 @@
 from __future__ import annotations
 
 from collections.abc import Sequence, Mapping
-from typing import Any
 from dataclasses import dataclass
-
-from qiskit.providers import Options
+from typing import Any
 
 
 @dataclass(frozen=True)
@@ -33,5 +31,5 @@ class StateFidelityResult:
     depending on the error mitigation method used."""
     metadata: Sequence[Mapping[str, Any]]
     """Additional information about the fidelity calculation."""
-    options: Options
-    """Primitive runtime options for the execution of the fidelity job."""
+    shots: int | Sequence[int]
+    """Primitive number of shots options for the execution of the fidelity job."""
diff --git a/qiskit_algorithms/time_evolvers/pvqd/pvqd.py b/qiskit_algorithms/time_evolvers/pvqd/pvqd.py
index c04d0bc3..8e04882b 100644
--- a/qiskit_algorithms/time_evolvers/pvqd/pvqd.py
+++ b/qiskit_algorithms/time_evolvers/pvqd/pvqd.py
@@ -20,7 +20,7 @@
 
 from qiskit.circuit import Parameter, ParameterVector, QuantumCircuit
 from qiskit.circuit.library import PauliEvolutionGate
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 from qiskit.synthesis import EvolutionSynthesis, LieTrotter
 from qiskit_algorithms.utils import algorithm_globals
@@ -74,7 +74,7 @@ class PVQD(RealTimeEvolver):
 
             from qiskit_algorithms.state_fidelities import ComputeUncompute
             from qiskit_algorithms.time_evolvers import TimeEvolutionProblem, PVQD
-            from qiskit.primitives import Estimator, Sampler
+            from qiskit.primitives import StatevectorEstimator as Estimator, StatevectorSampler as Sampler
             from qiskit.circuit.library import EfficientSU2
             from qiskit.quantum_info import SparsePauliOp, Pauli
             from qiskit_algorithms.optimizers import L_BFGS_B
@@ -121,7 +121,7 @@ def __init__(
         fidelity: BaseStateFidelity,
         ansatz: QuantumCircuit,
         initial_parameters: np.ndarray,
-        estimator: BaseEstimator | None = None,
+        estimator: BaseEstimatorV2 | None = None,
         optimizer: Optimizer | Minimizer | None = None,
         num_timesteps: int | None = None,
         evolution: EvolutionSynthesis | None = None,
diff --git a/qiskit_algorithms/time_evolvers/pvqd/utils.py b/qiskit_algorithms/time_evolvers/pvqd/utils.py
index b571a792..840626a3 100644
--- a/qiskit_algorithms/time_evolvers/pvqd/utils.py
+++ b/qiskit_algorithms/time_evolvers/pvqd/utils.py
@@ -21,7 +21,7 @@
 from qiskit.circuit import QuantumCircuit, Parameter, ParameterExpression
 from qiskit.compiler import transpile
 from qiskit.exceptions import QiskitError
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 
 from qiskit_algorithms.gradients import ParamShiftSamplerGradient as ParamShift
@@ -75,7 +75,7 @@ def _is_gradient_supported(ansatz: QuantumCircuit) -> bool:
 def _get_observable_evaluator(
     ansatz: QuantumCircuit,
     observables: BaseOperator | list[BaseOperator],
-    estimator: BaseEstimator,
+    estimator: BaseEstimatorV2,
 ) -> Callable[[np.ndarray], float | list[float]]:
     """Get a callable to evaluate a (list of) observable(s) for given circuit parameters."""
 
@@ -91,18 +91,10 @@ def evaluate_observables(theta: np.ndarray) -> float | list[float]:
         Raises:
             AlgorithmError: If a primitive job fails.
         """
-        if isinstance(observables, list):
-            num_observables = len(observables)
-            obs = observables
-        else:
-            num_observables = 1
-            obs = [observables]
-        states = [ansatz] * num_observables
-        parameter_values = [theta] * num_observables
 
         try:
-            estimator_job = estimator.run(states, obs, parameter_values=parameter_values)
-            results = estimator_job.result().values
+            estimator_job = estimator.run([(ansatz, observables, theta)])
+            results = estimator_job.result()[0].data.evs
         except Exception as exc:
             raise AlgorithmError("The primitive job failed!") from exc
 
diff --git a/qiskit_algorithms/time_evolvers/trotterization/trotter_qrte.py b/qiskit_algorithms/time_evolvers/trotterization/trotter_qrte.py
index 893cd985..4d216975 100644
--- a/qiskit_algorithms/time_evolvers/trotterization/trotter_qrte.py
+++ b/qiskit_algorithms/time_evolvers/trotterization/trotter_qrte.py
@@ -14,14 +14,17 @@
 
 from __future__ import annotations
 
+from typing import Any
+
 from qiskit import QuantumCircuit
 
 from qiskit.circuit.library import PauliEvolutionGate
 from qiskit.circuit.parametertable import ParameterView
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 from qiskit.quantum_info import Pauli, SparsePauliOp
 from qiskit.synthesis import ProductFormula, LieTrotter
 
+from qiskit_algorithms.custom_types import Transpiler
 from qiskit_algorithms.time_evolvers.time_evolution_problem import TimeEvolutionProblem
 from qiskit_algorithms.time_evolvers.time_evolution_result import TimeEvolutionResult
 from qiskit_algorithms.time_evolvers.real_time_evolver import RealTimeEvolver
@@ -41,7 +44,7 @@ class TrotterQRTE(RealTimeEvolver):
             from qiskit.quantum_info import Pauli, SparsePauliOp
             from qiskit import QuantumCircuit
             from qiskit_algorithms import TrotterQRTE, TimeEvolutionProblem
-            from qiskit.primitives import Estimator
+            from qiskit.primitives import StatevectorEstimator as Estimator
 
             operator = SparsePauliOp([Pauli("X"), Pauli("Z")])
             initial_state = QuantumCircuit(1)
@@ -56,8 +59,10 @@ class TrotterQRTE(RealTimeEvolver):
     def __init__(
         self,
         product_formula: ProductFormula | None = None,
-        estimator: BaseEstimator | None = None,
+        estimator: BaseEstimatorV2 | None = None,
         num_timesteps: int = 1,
+        transpiler: Transpiler | None = None,
+        transpiler_options: dict[str, Any] | None = None,
         *,
         insert_barriers: bool = False,
     ) -> None:
@@ -73,6 +78,11 @@ def __init__(
                 ``TimeEvolutionProblem.aux_operators``.
             num_timesteps: The number of time-steps the full evolution time is divided into
                 (repetitions of ``product_formula``).
+            transpiler: An optional object with a `run` method allowing to transpile the circuits
+                that are produced within this algorithm. If set to `None`, these won't be
+                transpiled.
+            transpiler_options: A dictionary of options to be passed to the transpiler's `run`
+                method as keyword arguments.
             insert_barriers: If True, insert a barrier after the initial state and after each Trotter
                 step.
         """
@@ -81,6 +91,8 @@ def __init__(
         self.num_timesteps = num_timesteps
         self.estimator = estimator
         self._insert_barriers = insert_barriers
+        self._transpiler = transpiler
+        self._transpiler_options = transpiler_options if transpiler_options is not None else {}
 
     @property
     def product_formula(self) -> ProductFormula:
@@ -96,14 +108,14 @@ def product_formula(self, product_formula: ProductFormula | None):
         self._product_formula = product_formula
 
     @property
-    def estimator(self) -> BaseEstimator | None:
+    def estimator(self) -> BaseEstimatorV2 | None:
         """
         Returns an estimator.
         """
         return self._estimator
 
     @estimator.setter
-    def estimator(self, estimator: BaseEstimator) -> None:
+    def estimator(self, estimator: BaseEstimatorV2) -> None:
         """
         Sets an estimator.
         """
@@ -197,6 +209,10 @@ def evolve(self, evolution_problem: TimeEvolutionProblem) -> TimeEvolutionResult
 
         evolved_state = QuantumCircuit(initial_state.num_qubits)
         evolved_state.append(initial_state, evolved_state.qubits)
+
+        if self._transpiler is not None:
+            evolved_state = self._transpiler.run(evolved_state, **self._transpiler_options)
+
         if self._insert_barriers:
             evolved_state.barrier()
 
@@ -235,6 +251,10 @@ def evolve(self, evolution_problem: TimeEvolutionProblem) -> TimeEvolutionResult
                     synthesis=self.product_formula,
                 )
             evolved_state.append(single_step_evolution_gate, evolved_state.qubits)
+
+            if self._transpiler is not None:
+                evolved_state = self._transpiler.run(evolved_state, **self._transpiler_options)
+
             if self._insert_barriers:
                 evolved_state.barrier()
 
diff --git a/qiskit_algorithms/time_evolvers/variational/var_qite.py b/qiskit_algorithms/time_evolvers/variational/var_qite.py
index a6f9d388..e2c720ff 100644
--- a/qiskit_algorithms/time_evolvers/variational/var_qite.py
+++ b/qiskit_algorithms/time_evolvers/variational/var_qite.py
@@ -21,7 +21,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.circuit import Parameter
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 
 from .solvers.ode.forward_euler_solver import ForwardEulerSolver
 
@@ -42,7 +42,7 @@ class VarQITE(VarQTE, ImaginaryTimeEvolver):
         from qiskit_algorithms.time_evolvers.variational import ImaginaryMcLachlanPrinciple
         from qiskit.circuit.library import EfficientSU2
         from qiskit.quantum_info import SparsePauliOp, Pauli
-        from qiskit.primitives import Estimator
+        from qiskit.primitives import StatevectorEstimator as Estimator
 
         observable = SparsePauliOp.from_list(
             [
@@ -78,7 +78,7 @@ def __init__(
         ansatz: QuantumCircuit,
         initial_parameters: Mapping[Parameter, float] | Sequence[float],
         variational_principle: ImaginaryVariationalPrinciple | None = None,
-        estimator: BaseEstimator | None = None,
+        estimator: BaseEstimatorV2 | None = None,
         ode_solver: Type[OdeSolver] | str = ForwardEulerSolver,
         lse_solver: Callable[[np.ndarray, np.ndarray], np.ndarray] | None = None,
         num_timesteps: int | None = None,
diff --git a/qiskit_algorithms/time_evolvers/variational/var_qrte.py b/qiskit_algorithms/time_evolvers/variational/var_qrte.py
index 12cda9a2..0b467296 100644
--- a/qiskit_algorithms/time_evolvers/variational/var_qrte.py
+++ b/qiskit_algorithms/time_evolvers/variational/var_qrte.py
@@ -21,7 +21,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.circuit import Parameter
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 
 from .solvers.ode.forward_euler_solver import ForwardEulerSolver
 
@@ -43,7 +43,7 @@ class VarQRTE(VarQTE, RealTimeEvolver):
         from qiskit_algorithms.time_evolvers.variational import RealMcLachlanPrinciple
         from qiskit.quantum_info import SparsePauliOp
         from qiskit.quantum_info import SparsePauliOp, Pauli
-        from qiskit.primitives import Estimator
+        from qiskit.primitives import StatevectorEstimator as Estimator
 
         observable = SparsePauliOp.from_list(
             [
@@ -79,7 +79,7 @@ def __init__(
         ansatz: QuantumCircuit,
         initial_parameters: Mapping[Parameter, float] | Sequence[float],
         variational_principle: RealVariationalPrinciple | None = None,
-        estimator: BaseEstimator | None = None,
+        estimator: BaseEstimatorV2 | None = None,
         ode_solver: Type[OdeSolver] | str = ForwardEulerSolver,
         lse_solver: Callable[[np.ndarray, np.ndarray], np.ndarray] | None = None,
         num_timesteps: int | None = None,
diff --git a/qiskit_algorithms/time_evolvers/variational/var_qte.py b/qiskit_algorithms/time_evolvers/variational/var_qte.py
index bc1b2e36..88939233 100644
--- a/qiskit_algorithms/time_evolvers/variational/var_qte.py
+++ b/qiskit_algorithms/time_evolvers/variational/var_qte.py
@@ -22,7 +22,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.circuit import Parameter
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives import BaseEstimatorV2
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 
 from .solvers.ode.forward_euler_solver import ForwardEulerSolver
@@ -77,7 +77,7 @@ def __init__(
         ansatz: QuantumCircuit,
         initial_parameters: Mapping[Parameter, float] | Sequence[float],
         variational_principle: VariationalPrinciple,
-        estimator: BaseEstimator,
+        estimator: BaseEstimatorV2,
         ode_solver: Type[OdeSolver] | str = ForwardEulerSolver,
         lse_solver: Callable[[np.ndarray, np.ndarray], np.ndarray] | None = None,
         num_timesteps: int | None = None,
diff --git a/qiskit_algorithms/time_evolvers/variational/variational_principles/imaginary_mc_lachlan_principle.py b/qiskit_algorithms/time_evolvers/variational/variational_principles/imaginary_mc_lachlan_principle.py
index bfa29efb..0560c5bd 100644
--- a/qiskit_algorithms/time_evolvers/variational/variational_principles/imaginary_mc_lachlan_principle.py
+++ b/qiskit_algorithms/time_evolvers/variational/variational_principles/imaginary_mc_lachlan_principle.py
@@ -21,7 +21,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.circuit import Parameter
-from qiskit.primitives import Estimator
+from qiskit.primitives import StatevectorEstimator as Estimator
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 
 from .imaginary_variational_principle import ImaginaryVariationalPrinciple
diff --git a/qiskit_algorithms/time_evolvers/variational/variational_principles/real_mc_lachlan_principle.py b/qiskit_algorithms/time_evolvers/variational/variational_principles/real_mc_lachlan_principle.py
index 822e6cc9..9a39277c 100644
--- a/qiskit_algorithms/time_evolvers/variational/variational_principles/real_mc_lachlan_principle.py
+++ b/qiskit_algorithms/time_evolvers/variational/variational_principles/real_mc_lachlan_principle.py
@@ -22,7 +22,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.circuit import Parameter
-from qiskit.primitives import Estimator
+from qiskit.primitives import StatevectorEstimator as Estimator
 from qiskit.quantum_info import SparsePauliOp
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 
@@ -103,8 +103,8 @@ def evolution_gradient(
         """
 
         try:
-            estimator_job = self.gradient._estimator.run([ansatz], [hamiltonian], [param_values])
-            energy = estimator_job.result().values[0]
+            estimator_job = self.gradient._estimator.run([(ansatz, hamiltonian, param_values)])
+            energy = estimator_job.result()[0].data.evs
         except Exception as exc:
             raise AlgorithmError("The primitive job failed!") from exc
 
diff --git a/test/eigensolvers/test_vqd.py b/test/eigensolvers/test_vqd.py
index c95d798c..33ad560f 100644
--- a/test/eigensolvers/test_vqd.py
+++ b/test/eigensolvers/test_vqd.py
@@ -20,7 +20,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.circuit.library import TwoLocal, RealAmplitudes
-from qiskit.primitives import Sampler, Estimator
+from qiskit.primitives import StatevectorSampler as Sampler, StatevectorEstimator as Estimator
 from qiskit.quantum_info import SparsePauliOp
 
 from qiskit_algorithms.eigensolvers import VQD, VQDResult
@@ -57,9 +57,8 @@ def setUp(self):
         )
         self.ry_wavefunction = TwoLocal(rotation_blocks="ry", entanglement_blocks="cz")
 
-        self.estimator = Estimator()
-        self.estimator_shots = Estimator(options={"shots": 1024, "seed": self.seed})
-        self.fidelity = ComputeUncompute(Sampler(options={"shots": 100_000, "seed": self.seed}))
+        self.estimator = Estimator(seed=self.seed)
+        self.fidelity = ComputeUncompute(Sampler(seed=self.seed, default_shots=10_000))
         self.betas = [3]
 
     @data(H2_SPARSE_PAULI)
@@ -92,11 +91,16 @@ def test_basic_operator(self, op):
 
         with self.subTest(msg="assert return ansatz is set"):
             job = self.estimator.run(
-                result.optimal_circuits,
-                [op] * len(result.optimal_points),
-                result.optimal_points,
+                [
+                    (circuits, op, optimal_points)
+                    for (circuits, optimal_points) in zip(
+                        result.optimal_circuits, result.optimal_points
+                    )
+                ]
             )
-            np.testing.assert_array_almost_equal(job.result().values, result.eigenvalues, 6)
+            job_result = job.result()
+            eigenvalues = np.array([job_result[i].data.evs for i in range(len(result.eigenvalues))])
+            np.testing.assert_array_almost_equal(eigenvalues, result.eigenvalues, 6)
 
         with self.subTest(msg="assert returned values are eigenvalues"):
             np.testing.assert_array_almost_equal(
@@ -123,9 +127,7 @@ def test_beta_autoeval(self, op):
         """Test beta auto-evaluation for different operator types."""
 
         with self.assertLogs(level="INFO") as logs:
-            vqd = VQD(
-                self.estimator_shots, self.fidelity, self.ryrz_wavefunction, optimizer=L_BFGS_B()
-            )
+            vqd = VQD(self.estimator, self.fidelity, self.ryrz_wavefunction, optimizer=L_BFGS_B())
             _ = vqd.compute_eigenvalues(op)
 
         # the first log message shows the value of beta[0]
@@ -179,7 +181,7 @@ def store_intermediate_result(eval_count, parameters, mean, metadata, step):
         wavefunction = self.ry_wavefunction
 
         vqd = VQD(
-            estimator=self.estimator_shots,
+            estimator=self.estimator,
             fidelity=self.fidelity,
             ansatz=wavefunction,
             optimizer=optimizer,
@@ -455,6 +457,18 @@ def test_convergence_threshold(self):
             SLSQP(),
             k=2,
             betas=self.betas,
+            initial_point=np.array(
+                [
+                    2.15707009,
+                    -2.6128808,
+                    1.40478697,
+                    -1.73909435,
+                    -2.89100903,
+                    1.75289926,
+                    -0.14760479,
+                    -2.00011645,
+                ]
+            ),
             convergence_threshold=1e-3,
         )
         with self.subTest("Failed convergence"):
@@ -465,7 +479,7 @@ def test_convergence_threshold(self):
             vqd.convergence_threshold = 1e-1
             result = vqd.compute_eigenvalues(operator=H2_SPARSE_PAULI)
             np.testing.assert_array_almost_equal(
-                result.eigenvalues.real, self.h2_energy_excited[:2], decimal=1
+                result.eigenvalues.real, self.h2_energy_excited[:2], decimal=2
             )
 
 
diff --git a/test/minimum_eigensolvers/test_adapt_vqe.py b/test/minimum_eigensolvers/test_adapt_vqe.py
index 4a13c9ca..011a3c93 100644
--- a/test/minimum_eigensolvers/test_adapt_vqe.py
+++ b/test/minimum_eigensolvers/test_adapt_vqe.py
@@ -21,7 +21,7 @@
 
 from qiskit.circuit import QuantumCircuit, QuantumRegister
 from qiskit.circuit.library import EvolvedOperatorAnsatz
-from qiskit.primitives import Estimator
+from qiskit.primitives import StatevectorEstimator
 from qiskit.quantum_info import SparsePauliOp
 
 from qiskit_algorithms.minimum_eigensolvers import VQE
@@ -83,7 +83,7 @@ def setUp(self):
 
     def test_default(self):
         """Default execution"""
-        calc = AdaptVQE(VQE(Estimator(), self.ansatz, self.optimizer))
+        calc = AdaptVQE(VQE(StatevectorEstimator(), self.ansatz, self.optimizer))
 
         res = calc.compute_minimum_eigenvalue(operator=self.h2_op)
 
@@ -98,7 +98,7 @@ def test_with_quantum_info(self):
             self.excitation_pool,
             initial_state=self.initial_state,
         )
-        calc = AdaptVQE(VQE(Estimator(), ansatz, self.optimizer))
+        calc = AdaptVQE(VQE(StatevectorEstimator(), ansatz, self.optimizer))
 
         res = calc.compute_minimum_eigenvalue(operator=self.h2_op)
 
@@ -110,7 +110,7 @@ def test_with_quantum_info(self):
     def test_converged(self):
         """Test to check termination criteria"""
         calc = AdaptVQE(
-            VQE(Estimator(), self.ansatz, self.optimizer),
+            VQE(StatevectorEstimator(), self.ansatz, self.optimizer),
             gradient_threshold=1e-3,
         )
         res = calc.compute_minimum_eigenvalue(operator=self.h2_op)
@@ -120,7 +120,7 @@ def test_converged(self):
     def test_maximum(self):
         """Test to check termination criteria"""
         calc = AdaptVQE(
-            VQE(Estimator(), self.ansatz, self.optimizer),
+            VQE(StatevectorEstimator(), self.ansatz, self.optimizer),
             max_iterations=1,
         )
         res = calc.compute_minimum_eigenvalue(operator=self.h2_op)
@@ -144,7 +144,7 @@ def test_eigenvalue_threshold(self):
         )
 
         calc = AdaptVQE(
-            VQE(Estimator(), ansatz, self.optimizer),
+            VQE(StatevectorEstimator(), ansatz, self.optimizer),
             eigenvalue_threshold=1,
         )
         res = calc.compute_minimum_eigenvalue(operator)
@@ -185,14 +185,14 @@ def test_cyclicity(self, seq, is_cycle):
 
     def test_vqe_solver(self):
         """Test to check if the VQE solver remains the same or not"""
-        solver = VQE(Estimator(), self.ansatz, self.optimizer)
+        solver = VQE(StatevectorEstimator(), self.ansatz, self.optimizer)
         calc = AdaptVQE(solver)
         _ = calc.compute_minimum_eigenvalue(operator=self.h2_op)
         self.assertEqual(solver.ansatz, calc.solver.ansatz)
 
     def test_gradient_calculation(self):
         """Test to check if the gradient calculation"""
-        solver = VQE(Estimator(), QuantumCircuit(1), self.optimizer)
+        solver = VQE(StatevectorEstimator(), QuantumCircuit(1), self.optimizer)
         calc = AdaptVQE(solver)
         calc._excitation_pool = [SparsePauliOp("X")]
         res = calc._compute_gradients(operator=SparsePauliOp("Y"), theta=[])
@@ -201,7 +201,7 @@ def test_gradient_calculation(self):
 
     def test_supports_aux_operators(self):
         """Test that auxiliary operators are supported"""
-        calc = AdaptVQE(VQE(Estimator(), self.ansatz, self.optimizer))
+        calc = AdaptVQE(VQE(StatevectorEstimator(), self.ansatz, self.optimizer))
         res = calc.compute_minimum_eigenvalue(operator=self.h2_op, aux_operators=[self.h2_op])
 
         expected_eigenvalue = -1.85727503
diff --git a/test/minimum_eigensolvers/test_qaoa.py b/test/minimum_eigensolvers/test_qaoa.py
index 17a15773..503cfa3e 100644
--- a/test/minimum_eigensolvers/test_qaoa.py
+++ b/test/minimum_eigensolvers/test_qaoa.py
@@ -23,7 +23,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.circuit import Parameter
-from qiskit.primitives import Sampler
+from qiskit.primitives import StatevectorSampler
 from qiskit.quantum_info import Pauli, SparsePauliOp
 from qiskit.result import QuasiDistribution
 
@@ -67,7 +67,7 @@ def setUp(self):
         super().setUp()
         self.seed = 10598
         algorithm_globals.random_seed = self.seed
-        self.sampler = Sampler()
+        self.sampler = StatevectorSampler()
 
     @idata(
         [
diff --git a/test/minimum_eigensolvers/test_sampling_vqe.py b/test/minimum_eigensolvers/test_sampling_vqe.py
index 35b84a3a..581b1b7e 100644
--- a/test/minimum_eigensolvers/test_sampling_vqe.py
+++ b/test/minimum_eigensolvers/test_sampling_vqe.py
@@ -23,7 +23,7 @@
 
 from qiskit.circuit import ParameterVector, QuantumCircuit
 from qiskit.circuit.library import RealAmplitudes, TwoLocal
-from qiskit.primitives import Sampler
+from qiskit.primitives import StatevectorSampler
 from qiskit.quantum_info import Pauli, SparsePauliOp
 
 from qiskit_algorithms import AlgorithmError
@@ -79,7 +79,7 @@ def test_exact_sampler(self, op):
         initial_point = np.zeros(ansatz.num_parameters)
         initial_point[-ansatz.num_qubits :] = np.pi / 2
 
-        vqe = SamplingVQE(Sampler(), ansatz, optimizer, initial_point=initial_point)
+        vqe = SamplingVQE(StatevectorSampler(), ansatz, optimizer, initial_point=initial_point)
         result = vqe.compute_minimum_eigenvalue(operator=op)
 
         with self.subTest(msg="test eigenvalue"):
@@ -107,7 +107,7 @@ def test_invalid_initial_point(self, op):
         ansatz = RealAmplitudes(2, reps=1)
         initial_point = np.array([1])
 
-        vqe = SamplingVQE(Sampler(), ansatz, SLSQP(), initial_point=initial_point)
+        vqe = SamplingVQE(StatevectorSampler(), ansatz, SLSQP(), initial_point=initial_point)
 
         with self.assertRaises(ValueError):
             _ = vqe.compute_minimum_eigenvalue(operator=op)
@@ -116,7 +116,7 @@ def test_invalid_initial_point(self, op):
     def test_ansatz_resize(self, op):
         """Test the ansatz is properly resized if it's a blueprint circuit."""
         ansatz = RealAmplitudes(1, reps=1)
-        vqe = SamplingVQE(Sampler(), ansatz, SLSQP())
+        vqe = SamplingVQE(StatevectorSampler(), ansatz, SLSQP())
         result = vqe.compute_minimum_eigenvalue(operator=op)
         self.assertAlmostEqual(result.eigenvalue, self.optimal_value, places=5)
 
@@ -125,7 +125,7 @@ def test_invalid_ansatz_size(self, op):
         """Test an error is raised if the ansatz has the wrong number of qubits."""
         ansatz = QuantumCircuit(1)
         ansatz.compose(RealAmplitudes(1, reps=2))
-        vqe = SamplingVQE(Sampler(), ansatz, SLSQP())
+        vqe = SamplingVQE(StatevectorSampler(), ansatz, SLSQP())
 
         with self.assertRaises(AlgorithmError):
             _ = vqe.compute_minimum_eigenvalue(operator=op)
@@ -134,7 +134,7 @@ def test_invalid_ansatz_size(self, op):
     def test_missing_varform_params(self, op):
         """Test specifying a variational form with no parameters raises an error."""
         circuit = QuantumCircuit(op.num_qubits)
-        vqe = SamplingVQE(Sampler(), circuit, SLSQP())
+        vqe = SamplingVQE(StatevectorSampler(), circuit, SLSQP())
         with self.assertRaises(AlgorithmError):
             vqe.compute_minimum_eigenvalue(operator=op)
 
@@ -142,7 +142,7 @@ def test_missing_varform_params(self, op):
     def test_batch_evaluate_slsqp(self, op):
         """Test batching with SLSQP (as representative of SciPyOptimizer)."""
         optimizer = SLSQP(max_evals_grouped=10)
-        vqe = SamplingVQE(Sampler(), RealAmplitudes(), optimizer)
+        vqe = SamplingVQE(StatevectorSampler(), RealAmplitudes(), optimizer)
         result = vqe.compute_minimum_eigenvalue(operator=op)
         self.assertAlmostEqual(result.eigenvalue, self.optimal_value, places=5)
 
@@ -150,8 +150,8 @@ def test_batch_evaluate_with_qnspsa(self):
         """Test batch evaluating with QNSPSA works."""
         ansatz = TwoLocal(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz")
 
-        wrapped_sampler = Sampler()
-        inner_sampler = Sampler()
+        wrapped_sampler = StatevectorSampler()
+        inner_sampler = StatevectorSampler()
 
         callcount = {"count": 0}
 
@@ -164,8 +164,7 @@ def wrapped_run(*args, **kwargs):
         fidelity = ComputeUncompute(wrapped_sampler)
 
         def fidelity_callable(left, right):
-            batchsize = np.asarray(left).shape[0]
-            job = fidelity.run(batchsize * [ansatz], batchsize * [ansatz], left, right)
+            job = fidelity.run(ansatz, ansatz, left, right)
             return job.result().fidelities
 
         qnspsa = QNSPSA(fidelity_callable, maxiter=5)
@@ -183,7 +182,7 @@ def fidelity_callable(left, right):
     def test_optimizer_scipy_callable(self):
         """Test passing a SciPy optimizer directly as callable."""
         vqe = SamplingVQE(
-            Sampler(),
+            StatevectorSampler(),
             RealAmplitudes(),
             partial(scipy_minimize, method="COBYLA", options={"maxiter": 2}),
         )
@@ -193,7 +192,7 @@ def test_optimizer_scipy_callable(self):
     def test_optimizer_callable(self):
         """Test passing a optimizer directly as callable."""
         ansatz = RealAmplitudes(1, reps=1)
-        vqe = SamplingVQE(Sampler(), ansatz, _mock_optimizer)
+        vqe = SamplingVQE(StatevectorSampler(), ansatz, _mock_optimizer)
         result = vqe.compute_minimum_eigenvalue(Pauli("Z"))
         self.assertTrue(np.all(result.optimal_point == np.zeros(ansatz.num_parameters)))
 
@@ -201,7 +200,7 @@ def test_optimizer_callable(self):
     def test_auxops(self, op):
         """Test passing auxiliary operators."""
         ansatz = RealAmplitudes(2, reps=1)
-        vqe = SamplingVQE(Sampler(), ansatz, SLSQP())
+        vqe = SamplingVQE(StatevectorSampler(), ansatz, SLSQP())
 
         as_list = [Pauli("ZZ"), Pauli("II")]
         with self.subTest(auxops=as_list):
@@ -220,7 +219,7 @@ def test_auxops(self, op):
 
     def test_nondiag_observable_raises(self):
         """Test passing a non-diagonal observable raises an error."""
-        vqe = SamplingVQE(Sampler(), RealAmplitudes(), SLSQP())
+        vqe = SamplingVQE(StatevectorSampler(), RealAmplitudes(), SLSQP())
 
         with self.assertRaises(ValueError):
             _ = vqe.compute_minimum_eigenvalue(Pauli("X"))
@@ -242,7 +241,7 @@ def store_intermediate_result(eval_count, parameters, mean, metadata):
             history["metadata"].append(metadata)
 
         sampling_vqe = SamplingVQE(
-            Sampler(),
+            StatevectorSampler(),
             RealAmplitudes(2, reps=1),
             SLSQP(),
             callback=store_intermediate_result,
@@ -271,7 +270,9 @@ def best_measurement(measurements):
 
         for aggregation in [alpha, best_measurement]:
             with self.subTest(aggregation=aggregation):
-                vqe = SamplingVQE(Sampler(), ansatz, _mock_optimizer, aggregation=best_measurement)
+                vqe = SamplingVQE(
+                    StatevectorSampler(), ansatz, _mock_optimizer, aggregation=best_measurement
+                )
                 result = vqe.compute_minimum_eigenvalue(Pauli("Z"))
 
                 # evaluation at x0=0 samples -1 and 1 with 50% probability, and our aggregation
diff --git a/test/minimum_eigensolvers/test_vqe.py b/test/minimum_eigensolvers/test_vqe.py
index 053cfd05..a1c64a8e 100644
--- a/test/minimum_eigensolvers/test_vqe.py
+++ b/test/minimum_eigensolvers/test_vqe.py
@@ -23,7 +23,7 @@
 from qiskit import QuantumCircuit
 from qiskit.circuit.library import RealAmplitudes, TwoLocal
 from qiskit.quantum_info import SparsePauliOp, Operator, Pauli
-from qiskit.primitives import Estimator, Sampler
+from qiskit.primitives import StatevectorEstimator, StatevectorSampler
 
 from qiskit_algorithms import AlgorithmError
 from qiskit_algorithms.gradients import ParamShiftEstimatorGradient
@@ -83,7 +83,7 @@ def setUp(self):
     @data(L_BFGS_B(), COBYLA())
     def test_using_ref_estimator(self, optimizer):
         """Test VQE using reference Estimator."""
-        vqe = VQE(Estimator(), self.ryrz_wavefunction, optimizer)
+        vqe = VQE(StatevectorEstimator(), self.ryrz_wavefunction, optimizer)
 
         result = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
 
@@ -109,9 +109,9 @@ def test_using_ref_estimator(self, optimizer):
             self.assertAlmostEqual(result.optimizer_result.fun, self.h2_energy, places=5)
 
         with self.subTest(msg="assert return ansatz is set"):
-            estimator = Estimator()
-            job = estimator.run(result.optimal_circuit, self.h2_op, result.optimal_point)
-            np.testing.assert_array_almost_equal(job.result().values, result.eigenvalue, 6)
+            estimator = StatevectorEstimator()
+            job = estimator.run([(result.optimal_circuit, self.h2_op, result.optimal_point)])
+            np.testing.assert_array_almost_equal(job.result()[0].data.evs, result.eigenvalue, 6)
 
     def test_invalid_initial_point(self):
         """Test the proper error is raised when the initial point has the wrong size."""
@@ -119,7 +119,7 @@ def test_invalid_initial_point(self):
         initial_point = np.array([1])
 
         vqe = VQE(
-            Estimator(),
+            StatevectorEstimator(),
             ansatz,
             SLSQP(),
             initial_point=initial_point,
@@ -131,7 +131,7 @@ def test_invalid_initial_point(self):
     def test_ansatz_resize(self):
         """Test the ansatz is properly resized if it's a blueprint circuit."""
         ansatz = RealAmplitudes(1, reps=1)
-        vqe = VQE(Estimator(), ansatz, SLSQP())
+        vqe = VQE(StatevectorEstimator(), ansatz, SLSQP())
         result = vqe.compute_minimum_eigenvalue(self.h2_op)
         self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
 
@@ -139,7 +139,7 @@ def test_invalid_ansatz_size(self):
         """Test an error is raised if the ansatz has the wrong number of qubits."""
         ansatz = QuantumCircuit(1)
         ansatz.compose(RealAmplitudes(1, reps=2))
-        vqe = VQE(Estimator(), ansatz, SLSQP())
+        vqe = VQE(StatevectorEstimator(), ansatz, SLSQP())
 
         with self.assertRaises(AlgorithmError):
             _ = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
@@ -147,7 +147,7 @@ def test_invalid_ansatz_size(self):
     def test_missing_ansatz_params(self):
         """Test specifying an ansatz with no parameters raises an error."""
         ansatz = QuantumCircuit(self.h2_op.num_qubits)
-        vqe = VQE(Estimator(), ansatz, SLSQP())
+        vqe = VQE(StatevectorEstimator(), ansatz, SLSQP())
         with self.assertRaises(AlgorithmError):
             vqe.compute_minimum_eigenvalue(operator=self.h2_op)
 
@@ -155,7 +155,7 @@ def test_max_evals_grouped(self):
         """Test with SLSQP with max_evals_grouped."""
         optimizer = SLSQP(maxiter=50, max_evals_grouped=5)
         vqe = VQE(
-            Estimator(),
+            StatevectorEstimator(),
             self.ryrz_wavefunction,
             optimizer,
         )
@@ -171,7 +171,7 @@ def test_max_evals_grouped(self):
     )
     def test_with_gradient(self, optimizer):
         """Test VQE using gradient primitive."""
-        estimator = Estimator()
+        estimator = StatevectorEstimator()
         vqe = VQE(
             estimator,
             self.ry_wavefunction,
@@ -184,7 +184,7 @@ def test_with_gradient(self, optimizer):
     def test_gradient_passed(self):
         """Test the gradient is properly passed into the optimizer."""
         inputs = {}
-        estimator = Estimator()
+        estimator = StatevectorEstimator()
         vqe = VQE(
             estimator,
             RealAmplitudes(),
@@ -197,7 +197,7 @@ def test_gradient_passed(self):
 
     def test_gradient_run(self):
         """Test using the gradient to calculate the minimum."""
-        estimator = Estimator()
+        estimator = StatevectorEstimator()
         vqe = VQE(
             estimator,
             RealAmplitudes(),
@@ -220,7 +220,7 @@ def store_intermediate_result(eval_count, parameters, mean, metadata):
         optimizer = COBYLA(maxiter=3)
         wavefunction = self.ry_wavefunction
 
-        estimator = Estimator()
+        estimator = StatevectorEstimator()
 
         vqe = VQE(
             estimator,
@@ -239,7 +239,7 @@ def store_intermediate_result(eval_count, parameters, mean, metadata):
     def test_reuse(self):
         """Test re-using a VQE algorithm instance."""
         ansatz = TwoLocal(rotation_blocks=["ry", "rz"], entanglement_blocks="cz")
-        vqe = VQE(Estimator(), ansatz, SLSQP(maxiter=300))
+        vqe = VQE(StatevectorEstimator(), ansatz, SLSQP(maxiter=300))
         with self.subTest(msg="assert VQE works once all info is available"):
             result = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
             self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
@@ -254,7 +254,7 @@ def test_reuse(self):
     def test_vqe_optimizer_reuse(self):
         """Test running same VQE twice to re-use optimizer, then switch optimizer"""
         vqe = VQE(
-            Estimator(),
+            StatevectorEstimator(),
             self.ryrz_wavefunction,
             SLSQP(),
         )
@@ -276,8 +276,8 @@ def test_default_batch_evaluation_on_spsa(self):
         """Test the default batching works."""
         ansatz = TwoLocal(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz")
 
-        wrapped_estimator = Estimator()
-        inner_estimator = Estimator()
+        wrapped_estimator = StatevectorEstimator()
+        inner_estimator = StatevectorEstimator()
 
         callcount = {"estimator": 0}
 
@@ -305,11 +305,11 @@ def test_batch_evaluate_with_qnspsa(self):
         """Test batch evaluating with QNSPSA works."""
         ansatz = TwoLocal(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz")
 
-        wrapped_sampler = Sampler()
-        inner_sampler = Sampler()
+        wrapped_sampler = StatevectorSampler()
+        inner_sampler = StatevectorSampler()
 
-        wrapped_estimator = Estimator()
-        inner_estimator = Estimator()
+        wrapped_estimator = StatevectorEstimator()
+        inner_estimator = StatevectorEstimator()
 
         callcount = {"sampler": 0, "estimator": 0}
 
@@ -328,7 +328,7 @@ def wrapped_sampler_run(*args, **kwargs):
 
         def fidelity_callable(left, right):
             batchsize = np.asarray(left).shape[0]
-            job = fidelity.run(batchsize * [ansatz], batchsize * [ansatz], left, right)
+            job = fidelity.run(ansatz, ansatz, left, right)
             return job.result().fidelities
 
         qnspsa = QNSPSA(fidelity_callable, maxiter=5)
@@ -353,7 +353,7 @@ def fidelity_callable(left, right):
     def test_optimizer_scipy_callable(self):
         """Test passing a SciPy optimizer directly as callable."""
         vqe = VQE(
-            Estimator(),
+            StatevectorEstimator(),
             self.ryrz_wavefunction,
             partial(scipy_minimize, method="L-BFGS-B", options={"maxiter": 10}),
         )
@@ -363,13 +363,13 @@ def test_optimizer_scipy_callable(self):
     def test_optimizer_callable(self):
         """Test passing a optimizer directly as callable."""
         ansatz = RealAmplitudes(1, reps=1)
-        vqe = VQE(Estimator(), ansatz, _mock_optimizer)
+        vqe = VQE(StatevectorEstimator(), ansatz, _mock_optimizer)
         result = vqe.compute_minimum_eigenvalue(SparsePauliOp("Z"))
         self.assertTrue(np.all(result.optimal_point == np.zeros(ansatz.num_parameters)))
 
     def test_aux_operators_list(self):
         """Test list-based aux_operators."""
-        vqe = VQE(Estimator(), self.ry_wavefunction, SLSQP(maxiter=300))
+        vqe = VQE(StatevectorEstimator(), self.ry_wavefunction, SLSQP(maxiter=300))
 
         with self.subTest("Test with an empty list."):
             result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators=[])
@@ -408,7 +408,7 @@ def test_aux_operators_list(self):
 
     def test_aux_operators_dict(self):
         """Test dictionary compatibility of aux_operators"""
-        vqe = VQE(Estimator(), self.ry_wavefunction, SLSQP(maxiter=300))
+        vqe = VQE(StatevectorEstimator(), self.ry_wavefunction, SLSQP(maxiter=300))
 
         with self.subTest("Test with an empty dictionary."):
             result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators={})
diff --git a/test/optimizers/test_optimizers.py b/test/optimizers/test_optimizers.py
index c3b9e4cb..595b9d25 100644
--- a/test/optimizers/test_optimizers.py
+++ b/test/optimizers/test_optimizers.py
@@ -13,6 +13,9 @@
 """Test Optimizers"""
 
 import unittest
+
+from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
+
 from test import QiskitAlgorithmsTestCase
 
 from typing import Optional, List, Tuple
@@ -23,7 +26,7 @@
 from qiskit.circuit.library import RealAmplitudes
 from qiskit.exceptions import MissingOptionalLibraryError
 from qiskit.utils import optionals
-from qiskit.primitives import Sampler
+from qiskit.primitives import StatevectorSampler as Sampler
 
 from qiskit_algorithms.optimizers import (
     ADAM,
@@ -407,7 +410,12 @@ def steps():
     def test_qnspsa(self):
         """Test QN-SPSA optimizer is serializable."""
         ansatz = RealAmplitudes(1)
-        fidelity = QNSPSA.get_fidelity(ansatz, sampler=Sampler())
+        fidelity = QNSPSA.get_fidelity(
+            ansatz,
+            sampler=Sampler(seed=123),
+            transpiler=generate_preset_pass_manager(optimization_level=1, seed_transpiler=42),
+            transpiler_options={"callable": lambda x: x},
+        )
         options = {
             "fidelity": fidelity,
             "maxiter": 100,
diff --git a/test/optimizers/test_spsa.py b/test/optimizers/test_spsa.py
index dca054e4..c77b2546 100644
--- a/test/optimizers/test_spsa.py
+++ b/test/optimizers/test_spsa.py
@@ -18,7 +18,8 @@
 import numpy as np
 
 from qiskit.circuit.library import PauliTwoDesign
-from qiskit.primitives import Estimator, Sampler
+from qiskit.primitives import StatevectorEstimator as Estimator, StatevectorSampler as Sampler
+
 from qiskit.quantum_info import SparsePauliOp, Statevector
 
 from qiskit_algorithms.optimizers import SPSA, QNSPSA
@@ -57,7 +58,7 @@ def objective(x):
             settings["regularization"] = 0.01
             expected_nfev = settings["maxiter"] * 5 + 1
         elif method == "qnspsa":
-            settings["fidelity"] = QNSPSA.get_fidelity(circuit, sampler=Sampler())
+            settings["fidelity"] = QNSPSA.get_fidelity(circuit, sampler=Sampler(seed=123))
             settings["regularization"] = 0.001
             settings["learning_rate"] = 0.05
             settings["perturbation"] = 0.05
@@ -204,7 +205,7 @@ def test_qnspsa_fidelity_primitives(self):
         initial_point = np.random.random(ansatz.num_parameters)
 
         with self.subTest(msg="pass as kwarg"):
-            fidelity = QNSPSA.get_fidelity(ansatz, sampler=Sampler())
+            fidelity = QNSPSA.get_fidelity(ansatz, sampler=Sampler(seed=123))
             result = fidelity(initial_point, initial_point)
 
             self.assertAlmostEqual(result[0], 1)
@@ -212,21 +213,19 @@ def test_qnspsa_fidelity_primitives(self):
     def test_qnspsa_max_evals_grouped(self):
         """Test using max_evals_grouped with QNSPSA."""
         circuit = PauliTwoDesign(3, reps=1, seed=1)
-        num_parameters = circuit.num_parameters
 
         obs = SparsePauliOp("ZZI")  # Z^Z^I
-        estimator = Estimator(options={"seed": 12})
+        estimator = Estimator(seed=12)
 
         initial_point = np.array(
             [0.82311034, 0.02611798, 0.21077064, 0.61842177, 0.09828447, 0.62013131]
         )
 
         def objective(x):
-            x = np.reshape(x, (-1, num_parameters)).tolist()
-            n = len(x)
-            return estimator.run(n * [circuit], n * [obs], x).result().values.real
+            results = estimator.run([(circuit, obs, x)]).result()
+            return np.array([res.data.evs for res in results]).real.reshape(-1)
 
-        fidelity = QNSPSA.get_fidelity(circuit, sampler=Sampler())
+        fidelity = QNSPSA.get_fidelity(circuit, sampler=Sampler(seed=12, default_shots=10_000))
         optimizer = QNSPSA(fidelity)
         optimizer.maxiter = 1
         optimizer.learning_rate = 0.05
diff --git a/test/state_fidelities/test_compute_uncompute.py b/test/state_fidelities/test_compute_uncompute.py
index a8cfe8f3..c95fef6f 100644
--- a/test/state_fidelities/test_compute_uncompute.py
+++ b/test/state_fidelities/test_compute_uncompute.py
@@ -13,17 +13,23 @@
 """Tests for Fidelity."""
 
 import unittest
+from itertools import product
+
+from ddt import data, ddt, unpack, idata
+from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
+
 from test import QiskitAlgorithmsTestCase
 
 import numpy as np
 
 from qiskit.circuit import QuantumCircuit, ParameterVector
 from qiskit.circuit.library import RealAmplitudes
-from qiskit.primitives import Sampler
+from qiskit.primitives import StatevectorSampler as Sampler
 
 from qiskit_algorithms.state_fidelities import ComputeUncompute
 
 
+@ddt
 class TestComputeUncompute(QiskitAlgorithmsTestCase):
     """Test Compute-Uncompute Fidelity class"""
 
@@ -49,7 +55,7 @@ def setUp(self):
         rx_rotation.h(1)
 
         self._circuit = [rx_rotations, ry_rotations, plus, zero, rx_rotation]
-        self._sampler = Sampler()
+        self._sampler = Sampler(seed=123, default_shots=10_000)
         self._left_params = np.array([[0, 0], [np.pi / 2, 0], [0, np.pi / 2], [np.pi, np.pi]])
         self._right_params = np.array([[0, 0], [0, 0], [np.pi / 2, 0], [0, 0]])
 
@@ -80,7 +86,7 @@ def test_local(self):
             job = fidelity.run(self._circuit[2], self._circuit[3])
             result = job.result()
             fidelities.append(result.fidelities[0])
-        np.testing.assert_allclose(fidelities, np.array([0.25, 0.5]), atol=1e-16)
+        np.testing.assert_allclose(fidelities, np.array([0.25, 0.5]), atol=1e-2, rtol=1e-2)
 
     def test_4param_pairs(self):
         """test for fidelity with four pairs of parameters"""
@@ -90,7 +96,9 @@ def test_4param_pairs(self):
             [self._circuit[0]] * n, [self._circuit[1]] * n, self._left_params, self._right_params
         )
         results = job.result()
-        np.testing.assert_allclose(results.fidelities, np.array([1.0, 0.5, 0.25, 0.0]), atol=1e-16)
+        np.testing.assert_allclose(
+            results.fidelities, np.array([1.0, 0.5, 0.25, 0.0]), atol=1e-2, rtol=1e-2
+        )
 
     def test_symmetry(self):
         """test for fidelity with the same circuit"""
@@ -111,11 +119,11 @@ def test_no_params(self):
         fidelity = ComputeUncompute(self._sampler)
         job = fidelity.run([self._circuit[2]], [self._circuit[3]])
         results = job.result()
-        np.testing.assert_allclose(results.fidelities, np.array([0.25]), atol=1e-16)
+        np.testing.assert_allclose(results.fidelities, np.array([0.25]), atol=1e-2, rtol=1e-2)
 
         job = fidelity.run([self._circuit[2]], [self._circuit[3]], [], [])
         results = job.result()
-        np.testing.assert_allclose(results.fidelities, np.array([0.25]), atol=1e-16)
+        np.testing.assert_allclose(results.fidelities, np.array([0.25]), atol=1e-2, rtol=1e-2)
 
     def test_left_param(self):
         """test for fidelity with only left parameters"""
@@ -125,7 +133,9 @@ def test_left_param(self):
             [self._circuit[1]] * n, [self._circuit[3]] * n, values_1=self._left_params
         )
         results = job.result()
-        np.testing.assert_allclose(results.fidelities, np.array([1.0, 0.5, 0.5, 0.0]), atol=1e-16)
+        np.testing.assert_allclose(
+            results.fidelities, np.array([1.0, 0.5, 0.5, 0.0]), atol=1e-2, rtol=1e-2
+        )
 
     def test_right_param(self):
         """test for fidelity with only right parameters"""
@@ -135,7 +145,9 @@ def test_right_param(self):
             [self._circuit[3]] * n, [self._circuit[1]] * n, values_2=self._left_params
         )
         results = job.result()
-        np.testing.assert_allclose(results.fidelities, np.array([1.0, 0.5, 0.5, 0.0]), atol=1e-16)
+        np.testing.assert_allclose(
+            results.fidelities, np.array([1.0, 0.5, 0.5, 0.0]), atol=1e-2, rtol=1e-2
+        )
 
     def test_not_set_circuits(self):
         """test for fidelity with no circuits."""
@@ -173,7 +185,9 @@ def test_asymmetric_params(self):
             [self._circuit[0]] * n, [self._circuit[4]] * n, self._left_params, right_params
         )
         result = job.result()
-        np.testing.assert_allclose(result.fidelities, np.array([0.5, 0.25, 0.25, 0.0]), atol=1e-16)
+        np.testing.assert_allclose(
+            result.fidelities, np.array([0.5, 0.25, 0.25, 0.0]), atol=1e-2, rtol=1e-2
+        )
 
     def test_input_format(self):
         """test for different input format variations"""
@@ -217,48 +231,63 @@ def test_input_measurements(self):
         result = job.result()
         np.testing.assert_allclose(result.fidelities, np.array([1.0]))
 
-    def test_options(self):
-        """Test fidelity's run options"""
-        sampler_shots = Sampler(options={"shots": 1024})
+    def test_shots(self):
+        """Test fidelity's run shots setting"""
+        sampler_shots = Sampler(default_shots=1024)
 
         with self.subTest("sampler"):
             # Only options in sampler
             fidelity = ComputeUncompute(sampler_shots)
-            options = fidelity.options
+            shots = fidelity.shots
             job = fidelity.run(self._circuit[2], self._circuit[3])
             result = job.result()
-            self.assertEqual(options.__dict__, {"shots": 1024})
-            self.assertEqual(result.options.__dict__, {"shots": 1024})
+            self.assertEqual(shots, None)
+            self.assertEqual(result.shots, 1024)
 
         with self.subTest("fidelity init"):
             # Fidelity default options override sampler
             # options and add new fields
-            fidelity = ComputeUncompute(sampler_shots, options={"shots": 2048, "dummy": 100})
-            options = fidelity.options
+            fidelity = ComputeUncompute(sampler_shots, shots=2048)
+            shots = fidelity.shots
             job = fidelity.run(self._circuit[2], self._circuit[3])
             result = job.result()
-            self.assertEqual(options.__dict__, {"shots": 2048, "dummy": 100})
-            self.assertEqual(result.options.__dict__, {"shots": 2048, "dummy": 100})
+            self.assertEqual(shots, 2048)
+            self.assertEqual(result.shots, 2048)
 
         with self.subTest("fidelity update"):
             # Update fidelity options
-            fidelity = ComputeUncompute(sampler_shots, options={"shots": 2048, "dummy": 100})
-            fidelity.update_default_options(shots=100)
-            options = fidelity.options
+            fidelity = ComputeUncompute(sampler_shots, shots=2048)
+            fidelity.shots = 100
+            shots = fidelity.shots
             job = fidelity.run(self._circuit[2], self._circuit[3])
             result = job.result()
-            self.assertEqual(options.__dict__, {"shots": 100, "dummy": 100})
-            self.assertEqual(result.options.__dict__, {"shots": 100, "dummy": 100})
+            self.assertEqual(shots, 100)
+            self.assertEqual(result.shots, 100)
 
         with self.subTest("fidelity run"):
             # Run options override fidelity options
-            fidelity = ComputeUncompute(sampler_shots, options={"shots": 2048, "dummy": 100})
-            job = fidelity.run(self._circuit[2], self._circuit[3], shots=50, dummy=None)
-            options = fidelity.options
+            fidelity = ComputeUncompute(sampler_shots, shots=2048)
+            job = fidelity.run(self._circuit[2], self._circuit[3], shots=50)
+            shots = fidelity.shots
             result = job.result()
             # Only default + sampler options. Not run.
-            self.assertEqual(options.__dict__, {"shots": 2048, "dummy": 100})
-            self.assertEqual(result.options.__dict__, {"shots": 50, "dummy": None})
+            self.assertEqual(shots, 2048)
+            self.assertEqual(result.shots, 50)
+
+    def test_transpiler(self):
+        """Test that the transpiler is called"""
+        pass_manager = generate_preset_pass_manager(optimization_level=1, seed_transpiler=42)
+        counts = [0]
+
+        def callback(**kwargs):
+            counts[0] = kwargs["count"]
+
+        fidelity = ComputeUncompute(
+            Sampler(), transpiler=pass_manager, transpiler_options={"callback": callback}
+        )
+        fidelity._construct_circuits(QuantumCircuit(1), QuantumCircuit(1))
+
+        self.assertEqual(counts[0], 15)
 
 
 if __name__ == "__main__":
diff --git a/test/test_grover.py b/test/test_grover.py
index 966c69a6..62f94dc2 100644
--- a/test/test_grover.py
+++ b/test/test_grover.py
@@ -12,20 +12,20 @@
 
 """Test Grover's algorithm."""
 
-import itertools
 import unittest
-from test import QiskitAlgorithmsTestCase
+from itertools import product
 
 import numpy as np
-from ddt import data, ddt, idata, unpack
-
+from ddt import data, ddt
 from qiskit import QuantumCircuit
 from qiskit.circuit.library import GroverOperator, PhaseOracle
-from qiskit.primitives import Sampler
+from qiskit.primitives import StatevectorSampler as Sampler
 from qiskit.quantum_info import Operator, Statevector
+from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
 from qiskit.utils.optionals import HAS_TWEEDLEDUM
 
 from qiskit_algorithms import AmplificationProblem, Grover
+from test import QiskitAlgorithmsTestCase
 
 
 @ddt
@@ -72,9 +72,7 @@ def is_good_state(bitstr):
                 # same as ``bitstr in ['01', '11']``
                 return bitstr[1] == "1"
 
-        possible_states = [
-            "".join(list(map(str, item))) for item in itertools.product([0, 1], repeat=2)
-        ]
+        possible_states = ["".join(list(map(str, item))) for item in product([0, 1], repeat=2)]
 
         oracle = QuantumCircuit(2)
         problem = AmplificationProblem(oracle, is_good_state=is_good_state)
@@ -91,50 +89,44 @@ class TestGrover(QiskitAlgorithmsTestCase):
 
     def setUp(self):
         super().setUp()
-        self._sampler = Sampler()
-        self._sampler_with_shots = Sampler(options={"shots": 1024, "seed": 123})
+        self._sampler = Sampler(seed=123)
 
     @unittest.skipUnless(HAS_TWEEDLEDUM, "tweedledum required for this test")
-    @data("ideal", "shots")
-    def test_implicit_phase_oracle_is_good_state(self, use_sampler):
+    def test_implicit_phase_oracle_is_good_state(self):
         """Test implicit default for is_good_state with PhaseOracle."""
-        grover = self._prepare_grover(use_sampler)
+        grover = self._prepare_grover()
         oracle = PhaseOracle("x & y")
         problem = AmplificationProblem(oracle)
         result = grover.amplify(problem)
         self.assertEqual(result.top_measurement, "11")
 
-    @idata(itertools.product(["ideal", "shots"], [[1, 2, 3], None, 2]))
-    @unpack
-    def test_iterations_with_good_state(self, use_sampler, iterations):
+    @data([1, 2, 3], None, 2)
+    def test_iterations_with_good_state(self, iterations):
         """Test the algorithm with different iteration types and with good state"""
-        grover = self._prepare_grover(use_sampler, iterations)
+        grover = self._prepare_grover(iterations)
         problem = AmplificationProblem(Statevector.from_label("111"), is_good_state=["111"])
         result = grover.amplify(problem)
         self.assertEqual(result.top_measurement, "111")
 
-    @idata(itertools.product(["shots"], [[1, 2, 3], None, 2]))
-    @unpack
-    def test_iterations_with_good_state_sample_from_iterations(self, use_sampler, iterations):
+    @data([1, 2, 3], None, 2)
+    def test_iterations_with_good_state_sample_from_iterations(self, iterations):
         """Test the algorithm with different iteration types and with good state"""
-        grover = self._prepare_grover(use_sampler, iterations, sample_from_iterations=True)
+        grover = self._prepare_grover(iterations, sample_from_iterations=True)
         problem = AmplificationProblem(Statevector.from_label("111"), is_good_state=["111"])
         result = grover.amplify(problem)
         self.assertEqual(result.top_measurement, "111")
 
-    @data("ideal", "shots")
-    def test_fixed_iterations_without_good_state(self, use_sampler):
+    def test_fixed_iterations_without_good_state(self):
         """Test the algorithm with iterations as an int and without good state"""
-        grover = self._prepare_grover(use_sampler, iterations=2)
+        grover = self._prepare_grover(iterations=2)
         problem = AmplificationProblem(Statevector.from_label("111"))
         result = grover.amplify(problem)
         self.assertEqual(result.top_measurement, "111")
 
-    @idata(itertools.product(["ideal", "shots"], [[1, 2, 3], None]))
-    @unpack
-    def test_iterations_without_good_state(self, use_sampler, iterations):
+    @data([1, 2, 3], None)
+    def test_iterations_without_good_state(self, iterations):
         """Test the correct error is thrown for none/list of iterations and without good state"""
-        grover = self._prepare_grover(use_sampler, iterations=iterations)
+        grover = self._prepare_grover(iterations=iterations)
         problem = AmplificationProblem(Statevector.from_label("111"))
 
         with self.assertRaisesRegex(
@@ -142,8 +134,7 @@ def test_iterations_without_good_state(self, use_sampler, iterations):
         ):
             grover.amplify(problem)
 
-    @data("ideal", "shots")
-    def test_iterator(self, use_sampler):
+    def test_iterator(self):
         """Test running the algorithm on an iterator."""
 
         # step-function iterator
@@ -155,63 +146,57 @@ def iterator():
                 if count % wait == 0:
                     value += 1
 
-        grover = self._prepare_grover(use_sampler, iterations=iterator())
+        grover = self._prepare_grover(iterations=iterator())
         problem = AmplificationProblem(Statevector.from_label("111"), is_good_state=["111"])
         result = grover.amplify(problem)
         self.assertEqual(result.top_measurement, "111")
 
-    @data("ideal", "shots")
-    def test_growth_rate(self, use_sampler):
+    def test_growth_rate(self):
         """Test running the algorithm on a growth rate"""
-        grover = self._prepare_grover(use_sampler, growth_rate=8 / 7)
+        grover = self._prepare_grover(growth_rate=8 / 7)
         problem = AmplificationProblem(Statevector.from_label("111"), is_good_state=["111"])
         result = grover.amplify(problem)
         self.assertEqual(result.top_measurement, "111")
 
-    @data("ideal", "shots")
-    def test_max_num_iterations(self, use_sampler):
+    def test_max_num_iterations(self):
         """Test the iteration stops when the maximum number of iterations is reached."""
 
         def zero():
             while True:
                 yield 0
 
-        grover = self._prepare_grover(use_sampler, iterations=zero())
+        grover = self._prepare_grover(iterations=zero())
         n = 5
         problem = AmplificationProblem(Statevector.from_label("1" * n), is_good_state=["1" * n])
         result = grover.amplify(problem)
         self.assertEqual(len(result.iterations), 2**n)
 
-    @data("ideal", "shots")
-    def test_max_power(self, use_sampler):
+    def test_max_power(self):
         """Test the iteration stops when the maximum power is reached."""
         lam = 10.0
-        grover = self._prepare_grover(use_sampler, growth_rate=lam)
+        grover = self._prepare_grover(growth_rate=lam)
         problem = AmplificationProblem(Statevector.from_label("111"), is_good_state=["111"])
         result = grover.amplify(problem)
         self.assertEqual(len(result.iterations), 0)
 
-    @data("ideal", "shots")
-    def test_run_circuit_oracle(self, use_sampler):
+    def test_run_circuit_oracle(self):
         """Test execution with a quantum circuit oracle"""
         oracle = QuantumCircuit(2)
         oracle.cz(0, 1)
         problem = AmplificationProblem(oracle, is_good_state=["11"])
-        grover = self._prepare_grover(use_sampler)
+        grover = self._prepare_grover()
         result = grover.amplify(problem)
         self.assertIn(result.top_measurement, ["11"])
 
-    @data("ideal", "shots")
-    def test_run_state_vector_oracle(self, use_sampler):
+    def test_run_state_vector_oracle(self):
         """Test execution with a state vector oracle"""
         mark_state = Statevector.from_label("11")
         problem = AmplificationProblem(mark_state, is_good_state=["11"])
-        grover = self._prepare_grover(use_sampler)
+        grover = self._prepare_grover()
         result = grover.amplify(problem)
         self.assertIn(result.top_measurement, ["11"])
 
-    @data("ideal", "shots")
-    def test_run_custom_grover_operator(self, use_sampler):
+    def test_run_custom_grover_operator(self):
         """Test execution with a grover operator oracle"""
         oracle = QuantumCircuit(2)
         oracle.cz(0, 1)
@@ -219,7 +204,7 @@ def test_run_custom_grover_operator(self, use_sampler):
         problem = AmplificationProblem(
             oracle=oracle, grover_operator=grover_op, is_good_state=["11"]
         )
-        grover = self._prepare_grover(use_sampler)
+        grover = self._prepare_grover()
         result = grover.amplify(problem)
         self.assertIn(result.top_measurement, ["11"])
 
@@ -247,14 +232,13 @@ def test_construct_circuit(self):
 
         self.assertTrue(Operator(constructed).equiv(Operator(expected)))
 
-    @data("ideal", "shots")
-    def test_circuit_result(self, use_sampler):
+    def test_circuit_result(self):
         """Test circuit_result"""
         oracle = QuantumCircuit(2)
         oracle.cz(0, 1)
         # is_good_state=['00'] is intentionally selected to obtain a list of results
         problem = AmplificationProblem(oracle, is_good_state=["00"])
-        grover = self._prepare_grover(use_sampler, iterations=[1, 2, 3, 4])
+        grover = self._prepare_grover(iterations=[1, 2, 3, 4])
 
         result = grover.amplify(problem)
 
@@ -267,23 +251,21 @@ def test_circuit_result(self, use_sampler):
                 self.assertTupleEqual(keys, ("00", "01", "10", "11"))
                 np.testing.assert_allclose(values, [0.25, 0.25, 0.25, 0.25], atol=0.2)
 
-    @data("ideal", "shots")
-    def test_max_probability(self, use_sampler):
+    def test_max_probability(self):
         """Test max_probability"""
         oracle = QuantumCircuit(2)
         oracle.cz(0, 1)
         problem = AmplificationProblem(oracle, is_good_state=["11"])
-        grover = self._prepare_grover(use_sampler)
+        grover = self._prepare_grover()
         result = grover.amplify(problem)
         self.assertAlmostEqual(result.max_probability, 1.0)
 
     @unittest.skipUnless(HAS_TWEEDLEDUM, "tweedledum required for this test")
-    @data("ideal", "shots")
-    def test_oracle_evaluation(self, use_sampler):
+    def test_oracle_evaluation(self):
         """Test oracle_evaluation for PhaseOracle"""
         oracle = PhaseOracle("x1 & x2 & (not x3)")
         problem = AmplificationProblem(oracle, is_good_state=oracle.evaluate_bitstring)
-        grover = self._prepare_grover(use_sampler)
+        grover = self._prepare_grover()
         result = grover.amplify(problem)
         self.assertTrue(result.oracle_evaluation)
         self.assertEqual("011", result.top_measurement)
@@ -294,27 +276,41 @@ def test_sampler_setter(self):
         grover.sampler = self._sampler
         self.assertEqual(grover.sampler, self._sampler)
 
+    def test_transpiler(self):
+        """Test that the transpiler is called"""
+        pass_manager = generate_preset_pass_manager(optimization_level=1, seed_transpiler=42)
+        counts = [0]
+
+        def callback(**kwargs):
+            counts[0] = kwargs["count"]
+
+        oracle = QuantumCircuit(2)
+        oracle.cz(0, 1)
+        # is_good_state=['00'] is intentionally selected to obtain a list of results
+        problem = AmplificationProblem(oracle)
+
+        Grover(
+            iterations=1,
+            sampler=Sampler(),
+            transpiler=pass_manager,
+            transpiler_options={"callback": callback},
+        ).amplify(problem)
+
+        self.assertEqual(counts[0], 15)
+
     def _prepare_grover(
-        self, use_sampler, iterations=None, growth_rate=None, sample_from_iterations=False
+        self,
+        iterations=None,
+        growth_rate=None,
+        sample_from_iterations=False,
     ):
         """Prepare Grover instance for test"""
-        if use_sampler == "ideal":
-            grover = Grover(
-                sampler=self._sampler,
-                iterations=iterations,
-                growth_rate=growth_rate,
-                sample_from_iterations=sample_from_iterations,
-            )
-        elif use_sampler == "shots":
-            grover = Grover(
-                sampler=self._sampler_with_shots,
-                iterations=iterations,
-                growth_rate=growth_rate,
-                sample_from_iterations=sample_from_iterations,
-            )
-        else:
-            raise RuntimeError("Unexpected `use_sampler` value {use_sampler}")
-        return grover
+        return Grover(
+            sampler=self._sampler,
+            iterations=iterations,
+            growth_rate=growth_rate,
+            sample_from_iterations=sample_from_iterations,
+        )
 
 
 if __name__ == "__main__":
diff --git a/test/test_phase_estimator.py b/test/test_phase_estimator.py
index d2edad1a..1f3e6e3b 100644
--- a/test/test_phase_estimator.py
+++ b/test/test_phase_estimator.py
@@ -11,23 +11,24 @@
 # that they have been altered from the originals.
 
 """Test phase estimation"""
-
 import unittest
-from test import QiskitAlgorithmsTestCase
-from ddt import ddt, data, unpack
+
 import numpy as np
-from qiskit.circuit.library import ZGate, XGate, HGate, IGate
-from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector, Operator
-from qiskit.synthesis import MatrixExponential, SuzukiTrotter
-from qiskit.primitives import Sampler
+from ddt import ddt, data, unpack
 from qiskit import QuantumCircuit
+from qiskit.circuit.library import HGate, XGate, IGate, ZGate
+from qiskit.primitives import StatevectorSampler as Sampler
+from qiskit.quantum_info import SparsePauliOp, Pauli, Statevector, Operator
+from qiskit.synthesis import MatrixExponential, SuzukiTrotter
+from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
 
-from qiskit_algorithms import PhaseEstimationScale
-from qiskit_algorithms.phase_estimators import (
-    PhaseEstimation,
+from qiskit_algorithms import (
     HamiltonianPhaseEstimation,
     IterativePhaseEstimation,
+    PhaseEstimation,
+    PhaseEstimationScale,
 )
+from test import QiskitAlgorithmsTestCase
 
 
 @ddt
@@ -45,9 +46,10 @@ def hamiltonian_pe_sampler(
         bound=None,
     ):
         """Run HamiltonianPhaseEstimation and return result with all phases."""
-        sampler = Sampler()
+        sampler = Sampler(default_shots=10_000, seed=42)
         phase_est = HamiltonianPhaseEstimation(
-            num_evaluation_qubits=num_evaluation_qubits, sampler=sampler
+            num_evaluation_qubits=num_evaluation_qubits,
+            sampler=sampler,
         )
         result = phase_est.estimate(
             hamiltonian=hamiltonian,
@@ -102,6 +104,21 @@ def test_single_pauli_op_sampler(self):
         with self.subTest("Second eigenvalue"):
             self.assertAlmostEqual(eigv, -0.98, delta=0.01)
 
+    def test_single_pauli_op_sampler_with_transpiler(self):
+        """Check that the transpilation does happen"""
+        pass_manager = generate_preset_pass_manager(optimization_level=1, seed_transpiler=42)
+        counts = [0]
+
+        def callback(**kwargs):
+            counts[0] = kwargs["count"]
+
+        phase_est = HamiltonianPhaseEstimation(
+            1, Sampler(), transpiler=pass_manager, transpiler_options={"callback": callback}
+        )
+        phase_est.estimate(hamiltonian=SparsePauliOp(Pauli("Z")))
+
+        self.assertEqual(counts[0], 15)
+
     @data(
         (Statevector(QuantumCircuit(2).compose(IGate()).compose(HGate()))),
         (QuantumCircuit(2).compose(IGate()).compose(HGate())),
@@ -164,11 +181,11 @@ def one_phase_sampler(
         """Run phase estimation with operator, eigenvalue pair `unitary_circuit`,
         `state_preparation`. Return the estimated phase as a value in :math:`[0,1)`.
         """
+
         if shots is not None:
-            options = {"shots": shots}
+            sampler = Sampler(default_shots=shots, seed=42)
         else:
-            options = {}
-        sampler = Sampler(options=options)
+            sampler = Sampler(seed=42)
         if phase_estimator is None:
             phase_estimator = IterativePhaseEstimation
         if phase_estimator == IterativePhaseEstimation:
@@ -211,11 +228,7 @@ def test_qpe_Z_sampler(self, state_preparation, expected_phase, shots, phase_est
     def test_qpe_X_plus_minus_sampler(self, state_preparation, expected_phase, phase_estimator):
         """eigenproblem X, (|+>, |->)"""
         unitary_circuit = QuantumCircuit(1).compose(XGate())
-        phase = self.one_phase_sampler(
-            unitary_circuit,
-            state_preparation,
-            phase_estimator,
-        )
+        phase = self.one_phase_sampler(unitary_circuit, state_preparation, phase_estimator)
         self.assertEqual(phase, expected_phase)
 
     @data(
@@ -230,11 +243,7 @@ def test_qpe_RZ_sampler(self, state_preparation, expected_phase, phase_estimator
         alpha = np.pi / 2
         unitary_circuit = QuantumCircuit(1)
         unitary_circuit.rz(alpha, 0)
-        phase = self.one_phase_sampler(
-            unitary_circuit,
-            state_preparation,
-            phase_estimator,
-        )
+        phase = self.one_phase_sampler(unitary_circuit, state_preparation, phase_estimator)
         self.assertEqual(phase, expected_phase)
 
     @data(
@@ -265,11 +274,7 @@ def test_qpe_two_qubit_unitary(self, state_preparation, expected_phase, phase_es
         unitary_circuit = QuantumCircuit(2)
         unitary_circuit.t(0)
         unitary_circuit.t(1)
-        phase = self.one_phase_sampler(
-            unitary_circuit,
-            state_preparation,
-            phase_estimator,
-        )
+        phase = self.one_phase_sampler(unitary_circuit, state_preparation, phase_estimator)
         self.assertEqual(phase, expected_phase)
 
     def test_check_num_iterations_sampler(self):
@@ -286,6 +291,36 @@ def test_phase_estimation_scale_from_operator(self):
         scale = PhaseEstimationScale.from_pauli_sum(op)
         self.assertEqual(scale._bound, 4.0)
 
+    @data((PhaseEstimation, 15), (IterativePhaseEstimation, 22))
+    @unpack
+    def test_transpiler(self, phase_estimator, expected_counts):
+        """Test that the transpiler is called"""
+        pass_manager = generate_preset_pass_manager(optimization_level=1, seed_transpiler=42)
+        counts = [0]
+
+        def callback(**kwargs):
+            counts[0] = kwargs["count"]
+
+        if phase_estimator == IterativePhaseEstimation:
+            p_est = IterativePhaseEstimation(
+                num_iterations=6,
+                sampler=Sampler(),
+                transpiler=pass_manager,
+                transpiler_options={"callback": callback},
+            )
+        elif phase_estimator == PhaseEstimation:
+            p_est = PhaseEstimation(
+                num_evaluation_qubits=6,
+                sampler=Sampler(),
+                transpiler=pass_manager,
+                transpiler_options={"callback": callback},
+            )
+        else:
+            raise ValueError("Unrecognized phase_estimator")
+
+        p_est.estimate(unitary=QuantumCircuit(1), state_preparation=None)
+        self.assertEqual(counts[0], expected_counts)
+
     # pylint: disable=too-many-positional-arguments
     def phase_estimation_sampler(
         self,
@@ -313,7 +348,7 @@ def test_qpe_Zplus_sampler(self, construct_circuit):
         """superposition eigenproblem Z, |+>"""
         unitary_circuit = QuantumCircuit(1).compose(ZGate())
         state_preparation = QuantumCircuit(1).compose(HGate())  # prepare |+>
-        sampler = Sampler()
+        sampler = Sampler(default_shots=10_000, seed=42)
         result = self.phase_estimation_sampler(
             unitary_circuit,
             sampler,
@@ -326,7 +361,7 @@ def test_qpe_Zplus_sampler(self, construct_circuit):
             self.assertEqual(list(phases.keys()), [0.0, 0.5])
 
         with self.subTest("test phases has correct probabilities"):
-            np.testing.assert_allclose(list(phases.values()), [0.5, 0.5])
+            np.testing.assert_allclose(list(phases.values()), [0.5, 0.5], atol=1e-2, rtol=1e-2)
 
         with self.subTest("test bitstring representation"):
             phases = result.filter_phases(1e-15, as_float=False)
diff --git a/test/time_evolvers/test_pvqd.py b/test/time_evolvers/test_pvqd.py
index 740ac69f..c8b02241 100644
--- a/test/time_evolvers/test_pvqd.py
+++ b/test/time_evolvers/test_pvqd.py
@@ -20,7 +20,7 @@
 
 from qiskit.circuit import Gate, Parameter, QuantumCircuit
 from qiskit.circuit.library import EfficientSU2
-from qiskit.primitives import Estimator, Sampler
+from qiskit.primitives import StatevectorEstimator as Estimator, StatevectorSampler as Sampler
 from qiskit.quantum_info import Pauli, SparsePauliOp
 
 from qiskit_algorithms import AlgorithmError
diff --git a/test/time_evolvers/test_trotter_qrte.py b/test/time_evolvers/test_trotter_qrte.py
index 3dda5710..9341420e 100644
--- a/test/time_evolvers/test_trotter_qrte.py
+++ b/test/time_evolvers/test_trotter_qrte.py
@@ -13,6 +13,9 @@
 """Test TrotterQRTE."""
 
 import unittest
+
+from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
+
 from test import QiskitAlgorithmsTestCase
 from ddt import ddt, data, unpack
 import numpy as np
@@ -23,7 +26,7 @@
 from qiskit.circuit.library import ZGate
 from qiskit.quantum_info import Statevector, Pauli, SparsePauliOp
 from qiskit.circuit import Parameter
-from qiskit.primitives import Estimator
+from qiskit.primitives import StatevectorEstimator as Estimator
 from qiskit.synthesis import SuzukiTrotter, QDrift
 
 from qiskit_algorithms import TimeEvolutionProblem, TrotterQRTE
@@ -239,6 +242,28 @@ def test_barriers(self, insert_barrier):
             expected_circuit.decompose(reps=3), evolution_result.evolved_state.decompose(reps=5)
         )
 
+    def test_transpiler(self):
+        """Test that the transpiler is called"""
+        pass_manager = generate_preset_pass_manager(optimization_level=1, seed_transpiler=42)
+        counts = [0]
+
+        def callback(**kwargs):
+            counts[0] = kwargs["count"]
+
+        operator = SparsePauliOp([Pauli("X"), Pauli("Z")])
+        initial_state = QuantumCircuit(1)
+        time = 1
+        evolution_problem = TimeEvolutionProblem(operator, time, initial_state)
+
+        trotter_qrte = TrotterQRTE(
+            estimator=Estimator(),
+            transpiler=pass_manager,
+            transpiler_options={"callback": callback},
+        )
+        trotter_qrte.evolve(evolution_problem)
+
+        self.assertEqual(counts[0], 15)
+
     # pylint: disable=too-many-positional-arguments
     @staticmethod
     def _run_error_test(initial_state, operator, aux_ops, estimator, t_param, param_value_dict):
diff --git a/test/time_evolvers/variational/test_var_qite.py b/test/time_evolvers/variational/test_var_qite.py
index 8ec64f54..f058b5eb 100644
--- a/test/time_evolvers/variational/test_var_qite.py
+++ b/test/time_evolvers/variational/test_var_qite.py
@@ -19,7 +19,7 @@
 
 from qiskit import QuantumCircuit
 from qiskit.circuit import Parameter
-from qiskit.primitives import Estimator
+from qiskit.primitives import StatevectorEstimator as Estimator
 from qiskit.quantum_info import SparsePauliOp, Pauli
 from qiskit.circuit.library import EfficientSU2
 from qiskit.quantum_info import Statevector
@@ -90,60 +90,28 @@ def test_run_d_1_with_aux_ops(self):
             1.939353810908912,
         ]
 
-        with self.subTest(msg="Test exact backend."):
-            algorithm_globals.random_seed = self.seed
-            estimator = Estimator()
-            qgt = LinCombQGT(estimator)
-            gradient = LinCombEstimatorGradient(estimator)
-            var_principle = ImaginaryMcLachlanPrinciple(qgt, gradient)
+        algorithm_globals.random_seed = self.seed
 
-            var_qite = VarQITE(
-                ansatz, init_param_values, var_principle, estimator, num_timesteps=25
-            )
-            evolution_result = var_qite.evolve(evolution_problem)
-
-            aux_ops = evolution_result.aux_ops_evaluated
-
-            parameter_values = evolution_result.parameter_values[-1]
-
-            expected_aux_ops = (-0.2177982985749799, 0.2556790598588627)
+        estimator = Estimator(seed=self.seed)
+        qgt = LinCombQGT(estimator)
+        gradient = LinCombEstimatorGradient(estimator)
+        var_principle = ImaginaryMcLachlanPrinciple(qgt, gradient)
 
-            for i, parameter_value in enumerate(parameter_values):
-                np.testing.assert_almost_equal(
-                    float(parameter_value), thetas_expected[i], decimal=2
-                )
+        var_qite = VarQITE(ansatz, init_param_values, var_principle, estimator, num_timesteps=25)
+        evolution_result = var_qite.evolve(evolution_problem)
 
-            np.testing.assert_array_almost_equal(
-                [result[0] for result in aux_ops], expected_aux_ops
-            )
+        aux_ops = evolution_result.aux_ops_evaluated
 
-        with self.subTest(msg="Test shot-based backend."):
-            algorithm_globals.random_seed = self.seed
+        parameter_values = evolution_result.parameter_values[-1]
 
-            estimator = Estimator(options={"shots": 4096, "seed": self.seed})
-            qgt = LinCombQGT(estimator)
-            gradient = LinCombEstimatorGradient(estimator)
-            var_principle = ImaginaryMcLachlanPrinciple(qgt, gradient)
+        expected_aux_ops = (-0.24629853310903974, 0.2518122871921184)
 
-            var_qite = VarQITE(
-                ansatz, init_param_values, var_principle, estimator, num_timesteps=25
+        for i, parameter_value in enumerate(parameter_values):
+            np.testing.assert_almost_equal(
+                float(parameter_value), thetas_expected_shots[i], decimal=2
             )
-            evolution_result = var_qite.evolve(evolution_problem)
-
-            aux_ops = evolution_result.aux_ops_evaluated
-
-            parameter_values = evolution_result.parameter_values[-1]
-
-            expected_aux_ops = (-0.24629853310903974, 0.2518122871921184)
-
-            for i, parameter_value in enumerate(parameter_values):
-                np.testing.assert_almost_equal(
-                    float(parameter_value), thetas_expected_shots[i], decimal=2
-                )
 
-            np.testing.assert_array_almost_equal(
-                [result[0] for result in aux_ops], expected_aux_ops
-            )
+        np.testing.assert_array_almost_equal([result[0] for result in aux_ops], expected_aux_ops)
 
     def test_run_d_1_t_7(self):
         """Test VarQITE for d = 1 and t = 7 with RK45 ODE solver."""
@@ -263,51 +231,27 @@ def test_run_d_1_time_dependent(self):
         state_expected = Statevector([0.34849948 + 0.0j, 0.93730897 + 0.0j]).to_dict()
         # the expected final state is Statevector([0.34849948+0.j, 0.93730897+0.j])
 
-        with self.subTest(msg="Test exact backend."):
-            algorithm_globals.random_seed = self.seed
-            estimator = Estimator()
-            var_principle = ImaginaryMcLachlanPrinciple()
-
-            var_qite = VarQITE(
-                ansatz, init_param_values, var_principle, estimator, num_timesteps=100
-            )
-            evolution_result = var_qite.evolve(evolution_problem)
-            evolved_state = evolution_result.evolved_state
-            parameter_values = evolution_result.parameter_values[-1]
-
-            for key, evolved_value in Statevector(evolved_state).to_dict().items():
-                # np.allclose works with complex numbers
-                self.assertTrue(np.allclose(evolved_value, state_expected[key], 1e-02))
+        algorithm_globals.random_seed = self.seed
 
-            for i, parameter_value in enumerate(parameter_values):
-                np.testing.assert_almost_equal(
-                    float(parameter_value), thetas_expected[i], decimal=2
-                )
-
-        with self.subTest(msg="Test shot-based backend."):
-            algorithm_globals.random_seed = self.seed
+        estimator = Estimator(seed=self.seed)
+        var_principle = ImaginaryMcLachlanPrinciple()
 
-            estimator = Estimator(options={"shots": 4 * 4096, "seed": self.seed})
-            var_principle = ImaginaryMcLachlanPrinciple()
+        var_qite = VarQITE(ansatz, init_param_values, var_principle, estimator, num_timesteps=100)
 
-            var_qite = VarQITE(
-                ansatz, init_param_values, var_principle, estimator, num_timesteps=100
-            )
-
-            evolution_result = var_qite.evolve(evolution_problem)
+        evolution_result = var_qite.evolve(evolution_problem)
 
-            evolved_state = evolution_result.evolved_state
+        evolved_state = evolution_result.evolved_state
 
-            parameter_values = evolution_result.parameter_values[-1]
+        parameter_values = evolution_result.parameter_values[-1]
 
-            for key, evolved_value in Statevector(evolved_state).to_dict().items():
-                # np.allclose works with complex numbers
-                self.assertTrue(np.allclose(evolved_value, state_expected[key], 1e-02))
+        for key, evolved_value in Statevector(evolved_state).to_dict().items():
+            # np.allclose works with complex numbers
+            self.assertTrue(np.allclose(evolved_value, state_expected[key], 1e-02))
 
-            for i, parameter_value in enumerate(parameter_values):
-                np.testing.assert_almost_equal(
-                    float(parameter_value), thetas_expected_shots[i], decimal=2
-                )
+        for i, parameter_value in enumerate(parameter_values):
+            np.testing.assert_almost_equal(
+                float(parameter_value), thetas_expected_shots[i], decimal=2
+            )
 
     # pylint: disable=too-many-positional-arguments
     def _test_helper(self, observable, thetas_expected, time, var_qite, decimal):
diff --git a/test/time_evolvers/variational/variational_principles/imaginary/test_imaginary_mc_lachlan_principle.py b/test/time_evolvers/variational/variational_principles/imaginary/test_imaginary_mc_lachlan_principle.py
index 9d8c54fa..b16460ed 100644
--- a/test/time_evolvers/variational/variational_principles/imaginary/test_imaginary_mc_lachlan_principle.py
+++ b/test/time_evolvers/variational/variational_principles/imaginary/test_imaginary_mc_lachlan_principle.py
@@ -23,7 +23,7 @@
 
 from qiskit.quantum_info import SparsePauliOp
 from qiskit.circuit.library import EfficientSU2
-from qiskit.primitives import Estimator
+from qiskit.primitives import StatevectorEstimator as Estimator
 
 from qiskit_algorithms.gradients import LinCombEstimatorGradient, DerivativeType
 from qiskit_algorithms.time_evolvers.variational import (
@@ -103,7 +103,7 @@ def test_calc_calc_evolution_gradient(self):
 
     def test_gradient_setting(self):
         """Test reactions to wrong gradient settings.."""
-        estimator = Estimator()
+        estimator = Estimator(seed=123)
         gradient = LinCombEstimatorGradient(estimator, derivative_type=DerivativeType.IMAG)
 
         with self.assertWarns(Warning):
diff --git a/test/time_evolvers/variational/variational_principles/real/test_real_mc_lachlan_principle.py b/test/time_evolvers/variational/variational_principles/real/test_real_mc_lachlan_principle.py
index 40195d16..60cef1ae 100644
--- a/test/time_evolvers/variational/variational_principles/real/test_real_mc_lachlan_principle.py
+++ b/test/time_evolvers/variational/variational_principles/real/test_real_mc_lachlan_principle.py
@@ -24,7 +24,7 @@
 
 from qiskit.quantum_info import SparsePauliOp
 from qiskit.circuit.library import EfficientSU2
-from qiskit.primitives import Estimator
+from qiskit.primitives import StatevectorEstimator as Estimator
 
 from qiskit_algorithms.gradients import LinCombEstimatorGradient, DerivativeType
 from qiskit_algorithms.time_evolvers.variational import (
@@ -108,7 +108,7 @@ def test_calc_evolution_gradient(self):
 
     def test_gradient_setting(self):
         """Test reactions to wrong gradient settings.."""
-        estimator = Estimator()
+        estimator = Estimator(seed=123)
         gradient = LinCombEstimatorGradient(estimator, derivative_type=DerivativeType.REAL)
 
         with self.assertWarns(Warning):