diff --git a/examples/python/ghz.py b/examples/python/ghz.py
index 469fdd268797..a92e913c7003 100644
--- a/examples/python/ghz.py
+++ b/examples/python/ghz.py
@@ -15,7 +15,7 @@
 """
 
 from qiskit import QuantumCircuit
-from qiskit import BasicAer, execute
+from qiskit import BasicAer, transpile
 
 
 ###############################################################
@@ -35,7 +35,7 @@
     qc.measure(i, i)
 
 sim_backend = BasicAer.get_backend("qasm_simulator")
-job = execute(qc, sim_backend, shots=1024)
+job = sim_backend.run(transpile(qc, sim_backend), shots=1024)
 result = job.result()
 print("Qasm simulator : ")
 print(result.get_counts(qc))
diff --git a/examples/python/initialize.py b/examples/python/initialize.py
index 945537c591be..30bb3c9c336d 100644
--- a/examples/python/initialize.py
+++ b/examples/python/initialize.py
@@ -15,7 +15,7 @@
 """
 
 import math
-from qiskit import QuantumCircuit, execute, BasicAer
+from qiskit import QuantumCircuit, transpile, BasicAer
 
 
 ###############################################################
@@ -59,7 +59,7 @@
 
 # Initialize on local simulator
 sim_backend = BasicAer.get_backend("qasm_simulator")
-job = execute(circuit, sim_backend, shots=shots)
+job = sim_backend.run(transpile(circuit, sim_backend), shots=shots)
 result = job.result()
 
 counts = result.get_counts(circuit)
diff --git a/examples/python/load_qasm.py b/examples/python/load_qasm.py
index d9da1f1dcd58..c4ef7d97422d 100644
--- a/examples/python/load_qasm.py
+++ b/examples/python/load_qasm.py
@@ -13,7 +13,7 @@
 """Example on how to load a file into a QuantumCircuit."""
 
 from qiskit import QuantumCircuit
-from qiskit import execute, BasicAer
+from qiskit import BasicAer
 
 circ = QuantumCircuit.from_qasm_file("examples/qasm/entangled_registers.qasm")
 print(circ)
@@ -23,7 +23,7 @@
 
 
 # Compile and run the Quantum circuit on a local simulator backend
-job_sim = execute(circ, sim_backend)
+job_sim = sim_backend.run(circ)
 sim_result = job_sim.result()
 
 # Show the results
diff --git a/examples/python/qft.py b/examples/python/qft.py
index 5849bb3813ae..467b8042b61f 100644
--- a/examples/python/qft.py
+++ b/examples/python/qft.py
@@ -16,7 +16,7 @@
 
 import math
 from qiskit import QuantumCircuit
-from qiskit import execute, BasicAer
+from qiskit import transpile, BasicAer
 
 
 ###############################################################
@@ -68,7 +68,7 @@ def qft(circ, n):
 
 print("Qasm simulator")
 sim_backend = BasicAer.get_backend("qasm_simulator")
-job = execute([qft3, qft4, qft5], sim_backend, shots=1024)
+job = sim_backend.run(transpile([qft3, qft4, qft5], sim_backend), shots=1024)
 result = job.result()
 print(result.get_counts(qft3))
 print(result.get_counts(qft4))
diff --git a/examples/python/rippleadd.py b/examples/python/rippleadd.py
index 72cec53dd171..6b02774b1c3c 100644
--- a/examples/python/rippleadd.py
+++ b/examples/python/rippleadd.py
@@ -17,7 +17,7 @@
 
 from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
 from qiskit import BasicAer
-from qiskit import execute
+from qiskit import transpile
 
 ###############################################################
 # Set the backend name and coupling map.
@@ -101,12 +101,12 @@ def unmajority(p, a, b, c):
 ###############################################################
 
 # First version: not mapped
-job = execute(qc, backend=backend, coupling_map=None, shots=1024)
+job = backend.run(transpile(qc, backend=backend, coupling_map=None), shots=1024)
 result = job.result()
 print(result.get_counts(qc))
 
 # Second version: mapped to 2x8 array coupling graph
-job = execute(qc, backend=backend, coupling_map=coupling_map, shots=1024)
+job = backend.run(transpile(qc, basis_gates=["u", "cx"], coupling_map=coupling_map), shots=1024)
 result = job.result()
 print(result.get_counts(qc))
 
diff --git a/examples/python/teleport.py b/examples/python/teleport.py
index 507a142ae9df..c7ca3559c4c7 100644
--- a/examples/python/teleport.py
+++ b/examples/python/teleport.py
@@ -19,7 +19,7 @@
 
 from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
 from qiskit import BasicAer
-from qiskit import execute
+from qiskit import transpile
 
 ###############################################################
 # Set the backend name and coupling map.
@@ -65,14 +65,17 @@
 
 # First version: not mapped
 initial_layout = {q[0]: 0, q[1]: 1, q[2]: 2}
-job = execute(qc, backend=backend, coupling_map=None, shots=1024, initial_layout=initial_layout)
+job = backend.run(
+    transpile(qc, backend=backend, coupling_map=None, initial_layout=initial_layout), shots=1024
+)
 
 result = job.result()
 print(result.get_counts(qc))
 
 # Second version: mapped to 2x8 array coupling graph
-job = execute(
-    qc, backend=backend, coupling_map=coupling_map, shots=1024, initial_layout=initial_layout
+job = backend.run(
+    transpile(qc, backend=backend, coupling_map=coupling_map, initial_layout=initial_layout),
+    shots=1024,
 )
 result = job.result()
 print(result.get_counts(qc))
diff --git a/examples/python/using_qiskit_terra_level_0.py b/examples/python/using_qiskit_terra_level_0.py
index ad20728341a5..ce46396c1d6d 100644
--- a/examples/python/using_qiskit_terra_level_0.py
+++ b/examples/python/using_qiskit_terra_level_0.py
@@ -11,10 +11,10 @@
 # that they have been altered from the originals.
 
 """
-Example showing how to use Qiskit-Terra at level 0 (novice).
+Example showing how to use Qiskit at level 0 (no optimization).
 
 This example shows the most basic way to user Terra. It builds some circuits
-and runs them on both the BasicAer (local Qiskit provider) or IBMQ (remote IBMQ provider).
+and runs them on both the BasicAer (local Qiskit provider) or IBM Quantum (remote IBM Quantum provider).
 
 To control the compile parameters we have provided a transpile function which can be used
 as a level 1 user.
@@ -23,7 +23,7 @@
 
 # Import the Qiskit modules
 from qiskit import QuantumCircuit
-from qiskit import execute, BasicAer
+from qiskit import transpile, BasicAer
 
 # making first circuit: bell state
 qc1 = QuantumCircuit(2, 2)
@@ -41,7 +41,8 @@
 print(BasicAer.backends())
 
 # running the job
-job_sim = execute([qc1, qc2], BasicAer.get_backend("qasm_simulator"))
+sim_backend = BasicAer.get_backend("qasm_simulator")
+job_sim = sim_backend.run(transpile([qc1, qc2], sim_backend))
 sim_result = job_sim.result()
 
 # Show the results
diff --git a/qiskit/execute_function.py b/qiskit/execute_function.py
index c50b6dcf2a50..b5aea06c0cd8 100644
--- a/qiskit/execute_function.py
+++ b/qiskit/execute_function.py
@@ -26,6 +26,7 @@
 from qiskit.providers.backend import Backend
 from qiskit.pulse import Schedule, ScheduleBlock
 from qiskit.exceptions import QiskitError
+from qiskit.utils.deprecation import deprecate_func
 
 logger = logging.getLogger(__name__)
 
@@ -35,6 +36,13 @@ def _log_submission_time(start_time, end_time):
     logger.info(log_msg)
 
 
+@deprecate_func(
+    additional_msg="This function combines ``transpile`` and ``backend.run``, which is covered "
+    "by ``Sampler`` :mod:`~qiskit.primitives`. Alternatively, you can also run :func:`.transpile` "
+    "followed by ``backend.run()``. The ``execute`` function is going to be removed "
+    "in 1.0.",
+    since="0.46.0",
+)
 def execute(
     experiments,
     backend,
diff --git a/qiskit/pulse/builder.py b/qiskit/pulse/builder.py
index 5ba542f24330..344604a41245 100644
--- a/qiskit/pulse/builder.py
+++ b/qiskit/pulse/builder.py
@@ -946,7 +946,7 @@ def build(
 
     .. code-block::
 
-        from qiskit import execute, pulse
+        from qiskit import transpile, pulse
         from qiskit.providers.fake_provider import FakeOpenPulse2Q
 
         backend = FakeOpenPulse2Q()
@@ -963,7 +963,7 @@ def build(
 
     .. code-block:: python
 
-        qiskit.execute(pulse_prog, backend)
+        backend.run(transpile(pulse_prog, backend))
 
     Args:
         backend (Backend): A Qiskit backend. If not supplied certain
diff --git a/qiskit/quantum_info/analysis/distance.py b/qiskit/quantum_info/analysis/distance.py
index 292ff77a706c..d24bb1f580df 100644
--- a/qiskit/quantum_info/analysis/distance.py
+++ b/qiskit/quantum_info/analysis/distance.py
@@ -76,7 +76,7 @@ def hellinger_fidelity(dist_p: dict, dist_q: dict) -> float:
 
         .. code-block::
 
-            from qiskit import QuantumCircuit, execute, BasicAer
+            from qiskit import QuantumCircuit, BasicAer
             from qiskit.quantum_info.analysis import hellinger_fidelity
 
             qc = QuantumCircuit(5, 5)
@@ -88,8 +88,8 @@ def hellinger_fidelity(dist_p: dict, dist_q: dict) -> float:
             qc.measure(range(5), range(5))
 
             sim = BasicAer.get_backend('qasm_simulator')
-            res1 = execute(qc, sim).result()
-            res2 = execute(qc, sim).result()
+            res1 = sim.run(qc).result()
+            res2 = sim.run(qc).result()
 
             hellinger_fidelity(res1.get_counts(), res2.get_counts())
 
diff --git a/qiskit/test/providers/backend.py b/qiskit/test/providers/backend.py
index 83832b8a29e7..6587c9ab167a 100644
--- a/qiskit/test/providers/backend.py
+++ b/qiskit/test/providers/backend.py
@@ -14,7 +14,7 @@
 
 from unittest import SkipTest
 
-from qiskit import execute
+from qiskit import transpile
 from ..base import QiskitTestCase
 from ..reference_circuits import ReferenceCircuits
 
@@ -69,7 +69,7 @@ def test_status(self):
 
     def test_run_circuit(self):
         """Test running a single circuit."""
-        job = execute(self.circuit, self.backend)
+        job = self.backend.run(transpile(self.circuit, self.backend))
         result = job.result()
         self.assertEqual(result.success, True)
         return result
diff --git a/releasenotes/notes/deprecate-execute-67a5c68b93722057.yaml b/releasenotes/notes/deprecate-execute-67a5c68b93722057.yaml
new file mode 100644
index 000000000000..5a0b56ea3b36
--- /dev/null
+++ b/releasenotes/notes/deprecate-execute-67a5c68b93722057.yaml
@@ -0,0 +1,30 @@
+---
+deprecations:
+  - |
+    Qiskit's :func:`~.execute_function.execute` function is deprecated.
+    This function served as a high-level wrapper around transpiling
+    a circuit with some transpile options and running it on a backend
+    with some run options. To do the same thing, you can explicitly
+    use the :func:`~.transpile` function (with appropriate transpile
+    options) followed by ``backend.run()`` (with appropriate run options).
+
+    For example, instead of running::
+
+      from qiskit import execute
+      job = execute(circuit, backend)
+
+    you can run::
+
+      from qiskit import transpile
+      new_circuit = transpile(circuit, backend)
+      job = backend.run(new_circuit)
+
+    Alternatively, the ``Sampler`` primitive is semantically equivalent to the
+    deprecated :func:`~.execute_function.execute` function. The class
+    :class:`.BackendSampler` is a generic wrapper for backends that do not support
+    primitives:
+
+      from qiskit.primitives import BackendSampler
+      sampler = BackendSampler(backend)
+      job = sampler.run(circuit)
+
diff --git a/test/python/algorithms/test_measure_error_mitigation.py b/test/python/algorithms/test_measure_error_mitigation.py
index 772f06007e13..6a92af1e233b 100644
--- a/test/python/algorithms/test_measure_error_mitigation.py
+++ b/test/python/algorithms/test_measure_error_mitigation.py
@@ -508,9 +508,9 @@ def test_tensor_subset_fitter(self):
         circuit.measure(2, 1)
         circuit.measure(0, 2)
 
-        result = execute(
-            circuit, backend, noise_model=noise_model, shots=1000, seed_simulator=0
-        ).result()
+        with self.assertWarns(DeprecationWarning):
+            job = execute(circuit, backend, noise_model=noise_model, shots=1000, seed_simulator=0)
+        result = job.result()
         with self.subTest(subset=subset):
             with self.assertWarns(DeprecationWarning):
                 new_result = fitter.subset_fitter([1, 2, 0]).filter.apply(result)
diff --git a/test/python/basicaer/test_basicaer_integration.py b/test/python/basicaer/test_basicaer_integration.py
index 0c698009350f..9dc00df843d1 100644
--- a/test/python/basicaer/test_basicaer_integration.py
+++ b/test/python/basicaer/test_basicaer_integration.py
@@ -16,7 +16,6 @@
 
 from qiskit import BasicAer
 from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
-from qiskit import execute
 from qiskit.result import Result
 from qiskit.providers.basicaer import BasicAerError
 from qiskit.test import QiskitTestCase
@@ -33,7 +32,7 @@ def setUp(self):
         self._qc2 = QuantumCircuit(qr, cr, name="qc2")
         self._qc1.measure(qr[0], cr[0])
         self.backend = BasicAer.get_backend("qasm_simulator")
-        self._result1 = execute(self._qc1, self.backend).result()
+        self._result1 = self.backend.run(self._qc1).result()
 
     def test_builtin_simulator_result_fields(self):
         """Test components of a result from a local simulator."""
@@ -52,7 +51,7 @@ def test_basicaer_execute(self):
         qc.cx(qubit_reg[0], qubit_reg[1])
         qc.measure(qubit_reg, clbit_reg)
 
-        job = execute(qc, self.backend)
+        job = self.backend.run(qc)
         result = job.result()
         self.assertIsInstance(result, Result)
 
@@ -66,7 +65,7 @@ def test_basicaer_execute_two(self):
         qc.measure(qubit_reg, clbit_reg)
         qc_extra = QuantumCircuit(qubit_reg, clbit_reg, name="extra")
         qc_extra.measure(qubit_reg, clbit_reg)
-        job = execute([qc, qc_extra], self.backend)
+        job = self.backend.run([qc, qc_extra])
         result = job.result()
         self.assertIsInstance(result, Result)
 
@@ -76,7 +75,7 @@ def test_basicaer_num_qubits(self):
         qc.x(0)
         qc.measure(0, 0)
         with self.assertRaises(BasicAerError):
-            execute(qc, self.backend)
+            self.backend.run(qc)
 
 
 if __name__ == "__main__":
diff --git a/test/python/basicaer/test_multi_registers_convention.py b/test/python/basicaer/test_multi_registers_convention.py
index ebefc22f5434..45abcd5453cf 100644
--- a/test/python/basicaer/test_multi_registers_convention.py
+++ b/test/python/basicaer/test_multi_registers_convention.py
@@ -12,7 +12,7 @@
 
 """Test executing multiple-register circuits on BasicAer."""
 
-from qiskit import BasicAer, execute
+from qiskit import BasicAer
 from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister
 from qiskit.quantum_info import Operator, Statevector, process_fidelity, state_fidelity
 from qiskit.test import QiskitTestCase
@@ -39,17 +39,17 @@ def test_circuit_multi(self):
 
         backend_sim = BasicAer.get_backend("qasm_simulator")
 
-        result = execute(qc, backend_sim, seed_transpiler=34342).result()
+        result = backend_sim.run(qc).result()
         counts = result.get_counts(qc)
 
         target = {"01 10": 1024}
 
         backend_sim = BasicAer.get_backend("statevector_simulator")
-        result = execute(circ, backend_sim, seed_transpiler=3438).result()
+        result = backend_sim.run(circ).result()
         state = result.get_statevector(circ)
 
         backend_sim = BasicAer.get_backend("unitary_simulator")
-        result = execute(circ, backend_sim, seed_transpiler=3438).result()
+        result = backend_sim.run(circ).result()
         unitary = Operator(result.get_unitary(circ))
 
         self.assertEqual(counts, target)
diff --git a/test/python/basicaer/test_qasm_simulator.py b/test/python/basicaer/test_qasm_simulator.py
index ca4d48d660fd..cb0f86761123 100644
--- a/test/python/basicaer/test_qasm_simulator.py
+++ b/test/python/basicaer/test_qasm_simulator.py
@@ -61,7 +61,8 @@ def assertExecuteLog(self, log_msg):
         qr = QuantumRegister(2, "qr")
         cr = ClassicalRegister(4, "cr")
         circuit = QuantumCircuit(qr, cr)
-        execute(circuit, backend=self.backend, shots=shots, seed_simulator=self.seed)
+        with self.assertWarns(DeprecationWarning):
+            execute(circuit, backend=self.backend, shots=shots, seed_simulator=self.seed)
         self.log_output.seek(0)
         # Filter unrelated log lines
         output_lines = self.log_output.readlines()
@@ -92,7 +93,7 @@ def test_measure_sampler_repeated_qubits(self):
         circuit.measure(qr[1], cr[2])
         circuit.measure(qr[0], cr[3])
         target = {"0110": shots}
-        job = execute(circuit, backend=self.backend, shots=shots, seed_simulator=self.seed)
+        job = self.backend.run(circuit, shots=shots, seed_simulator=self.seed)
         result = job.result()
         counts = result.get_counts(0)
         self.assertEqual(counts, target)
@@ -109,7 +110,7 @@ def test_measure_sampler_single_qubit(self):
             circuit.x(qr[qubit])
             circuit.measure(qr[qubit], cr[0])
             target = {"1": shots}
-            job = execute(circuit, backend=self.backend, shots=shots, seed_simulator=self.seed)
+            job = self.backend.run(circuit, shots=shots, seed_simulator=self.seed)
             result = job.result()
             counts = result.get_counts(0)
             self.assertEqual(counts, target)
@@ -146,7 +147,7 @@ def test_measure_sampler_partial_qubit(self):
         circuit.barrier(qr)
         circuit.measure(qr[3], cr[3])
         target = {"1011": shots}
-        job = execute(circuit, backend=self.backend, shots=shots, seed_simulator=self.seed)
+        job = self.backend.run(circuit, shots=shots, seed_simulator=self.seed)
         result = job.result()
         counts = result.get_counts(0)
         self.assertEqual(counts, target)
@@ -215,9 +216,8 @@ def test_if_statement(self):
         circuit_if_false.measure(qr[0], cr[0])
         circuit_if_false.measure(qr[1], cr[1])
         circuit_if_false.measure(qr[2], cr[2])
-        job = execute(
+        job = self.backend.run(
             [circuit_if_true, circuit_if_false],
-            backend=self.backend,
             shots=shots,
             seed_simulator=self.seed,
         )
@@ -252,7 +252,7 @@ def test_bit_cif_crossaffect(self):
         circuit.measure(qr[2], cr[2])
         circuit.h(qr[0]).c_if(cr[0], True)
         circuit.measure(qr[0], cr1[0])
-        job = execute(circuit, backend=self.backend, shots=shots, seed_simulator=self.seed)
+        job = self.backend.run(circuit, shots=shots, seed_simulator=self.seed)
         result = job.result().get_counts()
         target = {"0 110": 100}
         self.assertEqual(result, target)
@@ -292,7 +292,9 @@ def test_teleport(self):
         circuit.z(qr[2]).c_if(cr0, 1)
         circuit.x(qr[2]).c_if(cr1, 1)
         circuit.measure(qr[2], cr2[0])
-        job = execute(circuit, backend=self.backend, shots=shots, seed_simulator=self.seed)
+        job = self.backend.run(
+            transpile(circuit, self.backend), shots=shots, seed_simulator=self.seed
+        )
         results = job.result()
         data = results.get_counts("teleport")
         alice = {
@@ -346,7 +348,7 @@ def test_memory(self):
         circ.measure(qr[3], cr1[1])
 
         shots = 50
-        job = execute(circ, backend=self.backend, shots=shots, memory=True)
+        job = self.backend.run(circ, shots=shots, memory=True)
         result = job.result()
         memory = result.get_memory()
         self.assertEqual(len(memory), shots)
@@ -373,7 +375,7 @@ def test_unitary(self):
             circuit = QuantumCircuit(qr, cr)
             circuit.unitary(multi_x, qr)
             circuit.measure(qr, cr)
-            job = execute(circuit, self.backend, shots=shots)
+            job = self.backend.run(transpile(circuit), shots=shots)
             result = job.result()
             counts = result.get_counts(0)
             self.assertEqual(counts, target_counts)
diff --git a/test/python/basicaer/test_statevector_simulator.py b/test/python/basicaer/test_statevector_simulator.py
index 15561261cbf8..611fd86fd606 100644
--- a/test/python/basicaer/test_statevector_simulator.py
+++ b/test/python/basicaer/test_statevector_simulator.py
@@ -18,7 +18,7 @@
 from qiskit.providers.basicaer import StatevectorSimulatorPy
 from qiskit.test import ReferenceCircuits
 from qiskit.test import providers
-from qiskit import QuantumRegister, QuantumCircuit, execute
+from qiskit import QuantumRegister, QuantumCircuit, transpile
 from qiskit.quantum_info.random import random_unitary
 from qiskit.quantum_info import state_fidelity
 
@@ -76,7 +76,7 @@ def test_unitary(self):
                 # Simulate output on circuit
                 circuit = QuantumCircuit(qr)
                 circuit.unitary(unitary, qr)
-                job = execute(circuit, self.backend)
+                job = self.backend.run(transpile(circuit, self.backend))
                 result = job.result()
                 psi_out = result.get_statevector(0)
                 fidelity = state_fidelity(psi_target, psi_out)
diff --git a/test/python/basicaer/test_unitary_simulator.py b/test/python/basicaer/test_unitary_simulator.py
index 3930b7bbc247..13b504f2503f 100644
--- a/test/python/basicaer/test_unitary_simulator.py
+++ b/test/python/basicaer/test_unitary_simulator.py
@@ -16,7 +16,7 @@
 
 import numpy as np
 
-from qiskit import execute
+from qiskit import transpile
 from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister
 from qiskit.providers.basicaer import UnitarySimulatorPy
 from qiskit.quantum_info.operators.predicates import matrix_equal
@@ -35,7 +35,7 @@ class BasicAerUnitarySimulatorPyTest(providers.BackendTestCase):
     def test_basicaer_unitary_simulator_py(self):
         """Test unitary simulator."""
         circuits = self._test_circuits()
-        job = execute(circuits, backend=self.backend)
+        job = self.backend.run(transpile(circuits, self.backend))
         sim_unitaries = [job.result().get_unitary(circ) for circ in circuits]
         reference_unitaries = self._reference_unitaries()
         for u_sim, u_ref in zip(sim_unitaries, reference_unitaries):
@@ -108,7 +108,7 @@ def test_unitary(self):
                 # Simulate output on circuit
                 circuit = QuantumCircuit(qr)
                 circuit.unitary(unitary, qr)
-                job = execute(circuit, self.backend)
+                job = self.backend.run(transpile(circuit, self.backend))
                 result = job.result()
                 unitary_out = Operator(result.get_unitary(0))
                 fidelity = process_fidelity(unitary_target, unitary_out)
@@ -124,7 +124,7 @@ def test_global_phase(self):
         circuit.z(q[0])
         circuit.x(q[0])
 
-        job = execute(circuit, self.backend)
+        job = self.backend.run(transpile(circuit, self.backend))
         result = job.result()
 
         unitary_out = result.get_unitary(circuit)
diff --git a/test/python/circuit/library/test_functional_pauli_rotations.py b/test/python/circuit/library/test_functional_pauli_rotations.py
index 38f6568f117a..85bc16ec7a8e 100644
--- a/test/python/circuit/library/test_functional_pauli_rotations.py
+++ b/test/python/circuit/library/test_functional_pauli_rotations.py
@@ -18,7 +18,7 @@
 from ddt import ddt, data, unpack
 
 from qiskit.test.base import QiskitTestCase
-from qiskit import BasicAer, execute
+from qiskit import BasicAer, transpile
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import (
     LinearPauliRotations,
@@ -40,7 +40,7 @@ def assertFunctionIsCorrect(self, function_circuit, reference):
         circuit.append(function_circuit.to_instruction(), list(range(circuit.num_qubits)))
 
         backend = BasicAer.get_backend("statevector_simulator")
-        statevector = execute(circuit, backend).result().get_statevector()
+        statevector = backend.run(transpile(circuit, backend)).result().get_statevector()
 
         probabilities = defaultdict(float)
         for i, statevector_amplitude in enumerate(statevector):
diff --git a/test/python/circuit/library/test_integer_comparator.py b/test/python/circuit/library/test_integer_comparator.py
index 2c6459cf38e3..a98e3cca3cec 100644
--- a/test/python/circuit/library/test_integer_comparator.py
+++ b/test/python/circuit/library/test_integer_comparator.py
@@ -17,7 +17,7 @@
 from ddt import ddt, data, unpack
 
 from qiskit.test.base import QiskitTestCase
-from qiskit import BasicAer, execute
+from qiskit import BasicAer, transpile
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import IntegerComparator
 
@@ -34,7 +34,7 @@ def assertComparisonIsCorrect(self, comp, num_state_qubits, value, geq):
 
         # run simulation
         backend = BasicAer.get_backend("statevector_simulator")
-        statevector = execute(qc, backend).result().get_statevector()
+        statevector = backend.run(transpile(qc, backend)).result().get_statevector()
         for i, amplitude in enumerate(statevector):
             prob = np.abs(amplitude) ** 2
             if prob > 1e-6:
diff --git a/test/python/circuit/library/test_linear_amplitude_function.py b/test/python/circuit/library/test_linear_amplitude_function.py
index 5d00b35d41fd..b7c15bab2d9e 100644
--- a/test/python/circuit/library/test_linear_amplitude_function.py
+++ b/test/python/circuit/library/test_linear_amplitude_function.py
@@ -18,7 +18,7 @@
 from ddt import ddt, data, unpack
 import numpy as np
 
-from qiskit import BasicAer, execute
+from qiskit import BasicAer, transpile
 from qiskit.test.base import QiskitTestCase
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import LinearAmplitudeFunction
@@ -38,7 +38,7 @@ def assertFunctionIsCorrect(self, function_circuit, reference):
         circuit.append(function_circuit.to_instruction(), list(range(circuit.num_qubits)))
 
         backend = BasicAer.get_backend("statevector_simulator")
-        statevector = execute(circuit, backend).result().get_statevector()
+        statevector = backend.run(transpile(circuit, backend)).result().get_statevector()
 
         probabilities = defaultdict(float)
         for i, statevector_amplitude in enumerate(statevector):
diff --git a/test/python/circuit/library/test_phase_estimation.py b/test/python/circuit/library/test_phase_estimation.py
index a77b7ebec9ad..8a423128ad2a 100644
--- a/test/python/circuit/library/test_phase_estimation.py
+++ b/test/python/circuit/library/test_phase_estimation.py
@@ -17,7 +17,6 @@
 import numpy as np
 
 from qiskit.test.base import QiskitTestCase
-from qiskit import BasicAer, execute
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import PhaseEstimation, QFT
 from qiskit.quantum_info import Statevector
@@ -64,11 +63,7 @@ def assertPhaseEstimationIsCorrect(
             inplace=True,
         )
         circuit.compose(pec, inplace=True)
-        # TODO use Statevector for simulation once Qiskit/qiskit-terra#4681 is resolved
-        # actual = Statevector.from_instruction(circuit).data
-        backend = BasicAer.get_backend("statevector_simulator")
-        actual = execute(circuit, backend).result().get_statevector()
-
+        actual = Statevector.from_instruction(circuit).data
         np.testing.assert_almost_equal(reference, actual)
 
     def test_phase_estimation(self):
diff --git a/test/python/circuit/library/test_piecewise_chebyshev.py b/test/python/circuit/library/test_piecewise_chebyshev.py
index 4b8cb6054e1d..203f8448db53 100644
--- a/test/python/circuit/library/test_piecewise_chebyshev.py
+++ b/test/python/circuit/library/test_piecewise_chebyshev.py
@@ -18,7 +18,7 @@
 from ddt import ddt, data, unpack
 
 from qiskit.test.base import QiskitTestCase
-from qiskit import BasicAer, execute
+from qiskit import BasicAer, transpile
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library.arithmetic.piecewise_chebyshev import PiecewiseChebyshev
 
@@ -37,7 +37,7 @@ def assertFunctionIsCorrect(self, function_circuit, reference):
         circuit.append(function_circuit.to_instruction(), list(range(circuit.num_qubits)))
 
         backend = BasicAer.get_backend("statevector_simulator")
-        statevector = execute(circuit, backend).result().get_statevector()
+        statevector = backend.run(transpile(circuit, backend)).result().get_statevector()
 
         probabilities = defaultdict(float)
         for i, statevector_amplitude in enumerate(statevector):
diff --git a/test/python/circuit/library/test_weighted_adder.py b/test/python/circuit/library/test_weighted_adder.py
index 7bfe614e0bd6..c363222c2e00 100644
--- a/test/python/circuit/library/test_weighted_adder.py
+++ b/test/python/circuit/library/test_weighted_adder.py
@@ -18,7 +18,7 @@
 import numpy as np
 
 from qiskit.test.base import QiskitTestCase
-from qiskit import BasicAer, execute
+from qiskit import BasicAer, transpile
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import WeightedAdder
 
@@ -35,7 +35,7 @@ def assertSummationIsCorrect(self, adder):
         circuit.append(adder.to_instruction(), list(range(adder.num_qubits)))
 
         backend = BasicAer.get_backend("statevector_simulator")
-        statevector = execute(circuit, backend).result().get_statevector()
+        statevector = backend.run(transpile(circuit, backend)).result().get_statevector()
 
         probabilities = defaultdict(float)
         for i, statevector_amplitude in enumerate(statevector):
diff --git a/test/python/circuit/test_circuit_operations.py b/test/python/circuit/test_circuit_operations.py
index 3f773fb09abc..f7a3b0df5589 100644
--- a/test/python/circuit/test_circuit_operations.py
+++ b/test/python/circuit/test_circuit_operations.py
@@ -16,7 +16,7 @@
 import numpy as np
 from ddt import data, ddt
 
-from qiskit import BasicAer, ClassicalRegister, QuantumCircuit, QuantumRegister, execute
+from qiskit import BasicAer, ClassicalRegister, QuantumCircuit, QuantumRegister
 from qiskit.circuit import Gate, Instruction, Measure, Parameter, Barrier
 from qiskit.circuit.bit import Bit
 from qiskit.circuit.classicalregister import Clbit
@@ -188,7 +188,7 @@ def test_compose_circuit(self):
         qc3 = qc1.compose(qc2)
         backend = BasicAer.get_backend("qasm_simulator")
         shots = 1024
-        result = execute(qc3, backend=backend, shots=shots, seed_simulator=78).result()
+        result = backend.run(qc3, shots=shots, seed_simulator=78).result()
         counts = result.get_counts()
         target = {"00": shots / 2, "01": shots / 2}
         threshold = 0.04 * shots
@@ -210,7 +210,7 @@ def test_compose_circuit_and(self):
         qc3 = qc1 & qc2
         backend = BasicAer.get_backend("qasm_simulator")
         shots = 1024
-        result = execute(qc3, backend=backend, shots=shots, seed_simulator=78).result()
+        result = backend.run(qc3, shots=shots, seed_simulator=78).result()
         counts = result.get_counts()
         target = {"00": shots / 2, "01": shots / 2}
         threshold = 0.04 * shots
@@ -232,7 +232,7 @@ def test_compose_circuit_iand(self):
         qc1 &= qc2
         backend = BasicAer.get_backend("qasm_simulator")
         shots = 1024
-        result = execute(qc1, backend=backend, shots=shots, seed_simulator=78).result()
+        result = backend.run(qc1, shots=shots, seed_simulator=78).result()
         counts = result.get_counts()
         target = {"00": shots / 2, "01": shots / 2}
         threshold = 0.04 * shots
@@ -282,7 +282,7 @@ def test_tensor_circuit(self):
         qc3 = qc1.tensor(qc2)
         backend = BasicAer.get_backend("qasm_simulator")
         shots = 1024
-        result = execute(qc3, backend=backend, shots=shots, seed_simulator=78).result()
+        result = backend.run(qc3, shots=shots, seed_simulator=78).result()
         counts = result.get_counts()
         target = {"00": shots / 2, "01": shots / 2}
         threshold = 0.04 * shots
@@ -303,7 +303,7 @@ def test_tensor_circuit_xor(self):
         qc3 = qc1 ^ qc2
         backend = BasicAer.get_backend("qasm_simulator")
         shots = 1024
-        result = execute(qc3, backend=backend, shots=shots, seed_simulator=78).result()
+        result = backend.run(qc3, shots=shots, seed_simulator=78).result()
         counts = result.get_counts()
         target = {"00": shots / 2, "01": shots / 2}
         threshold = 0.04 * shots
@@ -324,7 +324,7 @@ def test_tensor_circuit_ixor(self):
         qc1 ^= qc2
         backend = BasicAer.get_backend("qasm_simulator")
         shots = 1024
-        result = execute(qc1, backend=backend, shots=shots, seed_simulator=78).result()
+        result = backend.run(qc1, shots=shots, seed_simulator=78).result()
         counts = result.get_counts()
         target = {"00": shots / 2, "01": shots / 2}
         threshold = 0.04 * shots
diff --git a/test/python/circuit/test_controlled_gate.py b/test/python/circuit/test_controlled_gate.py
index 34bf21f30e2e..39042ae3d3ac 100644
--- a/test/python/circuit/test_controlled_gate.py
+++ b/test/python/circuit/test_controlled_gate.py
@@ -19,7 +19,7 @@
 from numpy import pi
 from ddt import ddt, data, unpack
 
-from qiskit import QuantumRegister, QuantumCircuit, execute, BasicAer, QiskitError
+from qiskit import QuantumRegister, QuantumCircuit, BasicAer, QiskitError, transpile
 from qiskit.test import QiskitTestCase
 from qiskit.circuit import ControlledGate, Parameter, Gate
 from qiskit.circuit.singleton import SingletonControlledGate, _SingletonControlledGateOverrides
@@ -424,10 +424,11 @@ def test_multi_control_u1(self):
         c_cu1 = cu1gate.control(1)
         qc_cu1.append(c_cu1, qr, [])
 
-        job = execute(
-            [qcnu1, qu1, qcu1, qc_cu1],
-            BasicAer.get_backend("unitary_simulator"),
-            basis_gates=["u1", "u2", "u3", "id", "cx"],
+        backend = BasicAer.get_backend("unitary_simulator")
+        job = backend.run(
+            transpile(
+                [qcnu1, qu1, qcu1, qc_cu1], backend, basis_gates=["u1", "u2", "u3", "id", "cx"]
+            ),
         )
         result = job.result()
 
@@ -487,8 +488,7 @@ def test_multi_controlled_u1_matrix(self, num_controls):
                 if bit == "0":
                     qc.x(q_controls[idx])
 
-            backend = BasicAer.get_backend("unitary_simulator")
-            simulated = execute(qc, backend).result().get_unitary(qc)
+            simulated = Operator(qc)
 
             base = PhaseGate(lam).to_matrix()
             expected = _compute_control_matrix(base, num_controls, ctrl_state=ctrl_state)
@@ -520,7 +520,7 @@ def test_multi_control_toffoli_matrix_clean_ancillas(self, num_controls):
 
         # execute the circuit and obtain statevector result
         backend = BasicAer.get_backend("unitary_simulator")
-        simulated = execute(qc, backend).result().get_unitary(qc)
+        simulated = backend.run(transpile(qc, backend)).result().get_unitary(qc)
 
         # compare to expectation
         if num_ancillas > 0:
@@ -551,7 +551,8 @@ def test_multi_control_toffoli_matrix_basic_dirty_ancillas(self, num_controls):
 
         qc.mcx(q_controls, q_target[0], q_ancillas, mode="basic-dirty-ancilla")
 
-        simulated = execute(qc, BasicAer.get_backend("unitary_simulator")).result().get_unitary(qc)
+        simulator = BasicAer.get_backend("unitary_simulator")
+        simulated = simulator.run(transpile(qc, simulator)).result().get_unitary(qc)
         if num_ancillas > 0:
             simulated = simulated[: 2 ** (num_controls + 1), : 2 ** (num_controls + 1)]
 
@@ -580,7 +581,8 @@ def test_multi_control_toffoli_matrix_advanced_dirty_ancillas(self, num_controls
 
         qc.mcx(q_controls, q_target[0], q_ancillas, mode="advanced")
 
-        simulated = execute(qc, BasicAer.get_backend("unitary_simulator")).result().get_unitary(qc)
+        simulator = BasicAer.get_backend("unitary_simulator")
+        simulated = simulator.run(transpile(qc, simulator)).result().get_unitary(qc)
         if num_ancillas > 0:
             simulated = simulated[: 2 ** (num_controls + 1), : 2 ** (num_controls + 1)]
 
@@ -601,7 +603,8 @@ def test_multi_control_toffoli_matrix_noancilla_dirty_ancillas(self, num_control
 
         qc.mcx(q_controls, q_target[0], None, mode="noancilla")
 
-        simulated = execute(qc, BasicAer.get_backend("unitary_simulator")).result().get_unitary(qc)
+        simulator = BasicAer.get_backend("unitary_simulator")
+        simulated = simulator.run(transpile(qc, simulator)).result().get_unitary(qc)
 
         base = XGate().to_matrix()
         expected = _compute_control_matrix(base, num_controls)
@@ -663,7 +666,7 @@ def test_multi_controlled_rotation_gate_matrices(
                     self.assertTrue(gates_used.issubset({"x", "u", "p", "cx"}))
 
             backend = BasicAer.get_backend("unitary_simulator")
-            simulated = execute(qc, backend).result().get_unitary(qc)
+            simulated = backend.run(transpile(qc, backend)).result().get_unitary(qc)
 
             if base_gate_name == "x":
                 rot_mat = RXGate(theta).to_matrix()
@@ -723,7 +726,7 @@ def test_multi_controlled_y_rotation_matrix_basic_mode(self, num_controls, use_b
             rot_mat = RYGate(theta).to_matrix()
 
             backend = BasicAer.get_backend("unitary_simulator")
-            simulated = execute(qc, backend).result().get_unitary(qc)
+            simulated = backend.run(transpile(qc, backend)).result().get_unitary(qc)
             if num_ancillas > 0:
                 simulated = simulated[: 2 ** (num_controls + 1), : 2 ** (num_controls + 1)]
 
@@ -778,7 +781,7 @@ def test_mcx_gates(self, num_ctrl_qubits):
                 if num_ctrl_qubits > 0:
                     circuit.x(list(range(num_ctrl_qubits)))
                 circuit.append(gate, list(range(gate.num_qubits)), [])
-                statevector = execute(circuit, backend).result().get_statevector()
+                statevector = backend.run(transpile(circuit, backend)).result().get_statevector()
 
                 # account for ancillas
                 if hasattr(gate, "num_ancilla_qubits") and gate.num_ancilla_qubits > 0:
@@ -1081,7 +1084,7 @@ def test_relative_phase_toffoli_gates(self, num_ctrl_qubits):
         else:  # num_ctrl_qubits == 3:
             circuit.rcccx(0, 1, 2, 3)
         simulator = BasicAer.get_backend("unitary_simulator")
-        simulated_mat = execute(circuit, simulator).result().get_unitary()
+        simulated_mat = simulator.run(transpile(circuit, simulator)).result().get_unitary()
 
         # get the matrix representation from the class itself
         if num_ctrl_qubits == 2:
diff --git a/test/python/circuit/test_diagonal_gate.py b/test/python/circuit/test_diagonal_gate.py
index fcf2b832a72d..eb73a9605b67 100644
--- a/test/python/circuit/test_diagonal_gate.py
+++ b/test/python/circuit/test_diagonal_gate.py
@@ -16,7 +16,7 @@
 import unittest
 import numpy as np
 
-from qiskit import QuantumCircuit, QuantumRegister, BasicAer, execute, assemble
+from qiskit import QuantumCircuit, QuantumRegister, BasicAer, assemble
 
 from qiskit import QiskitError
 from qiskit.test import QiskitTestCase
@@ -52,7 +52,7 @@ def test_diag_gate(self):
                 qc = transpile(qc, basis_gates=["u1", "u3", "u2", "cx", "id"], optimization_level=0)
                 # Simulate the decomposed gate
                 simulator = BasicAer.get_backend("unitary_simulator")
-                result = execute(qc, simulator).result()
+                result = simulator.run(qc).result()
                 unitary = result.get_unitary(qc)
                 unitary_desired = _get_diag_gate_matrix(diag)
                 self.assertTrue(matrix_equal(unitary, unitary_desired, ignore_phase=False))
diff --git a/test/python/circuit/test_extensions_standard.py b/test/python/circuit/test_extensions_standard.py
index 9e43fccdf1da..9f69d13400d6 100644
--- a/test/python/circuit/test_extensions_standard.py
+++ b/test/python/circuit/test_extensions_standard.py
@@ -19,7 +19,7 @@
 import numpy as np
 from scipy.linalg import expm
 from ddt import data, ddt, unpack
-from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister, execute
+from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister
 from qiskit.exceptions import QiskitError
 from qiskit.circuit.exceptions import CircuitError
 from qiskit.test import QiskitTestCase
@@ -37,7 +37,6 @@
     YGate,
     GlobalPhaseGate,
 )
-from qiskit import BasicAer
 from qiskit.quantum_info import Pauli
 from qiskit.quantum_info.operators.predicates import matrix_equal, is_unitary_matrix
 from qiskit.utils.optionals import HAS_TWEEDLEDUM
@@ -1380,7 +1379,6 @@ def test_to_matrix(self):
 
         params = [0.1 * (i + 1) for i in range(10)]
         gate_class_list = Gate.__subclasses__() + ControlledGate.__subclasses__()
-        simulator = BasicAer.get_backend("unitary_simulator")
         for gate_class in gate_class_list:
             if hasattr(gate_class, "__abstractmethods__"):
                 # gate_class is abstract
@@ -1410,11 +1408,10 @@ def test_to_matrix(self):
                 # gate doesn't implement to_matrix method: skip
                 self.log.info('to_matrix method FAILED for "%s" gate', gate.name)
                 continue
-            definition_unitary = execute([circ], simulator).result().get_unitary()
+            definition_unitary = Operator(circ).data
 
             with self.subTest(gate_class):
-                # TODO check for exact equality once BasicAer can handle global phase
-                self.assertTrue(matrix_equal(definition_unitary, gate_matrix, ignore_phase=True))
+                self.assertTrue(matrix_equal(definition_unitary, gate_matrix))
                 self.assertTrue(is_unitary_matrix(gate_matrix))
 
     @unittest.skipUnless(HAS_TWEEDLEDUM, "tweedledum required for this test")
diff --git a/test/python/circuit/test_initializer.py b/test/python/circuit/test_initializer.py
index 4710cfcaeba0..23666297a3e6 100644
--- a/test/python/circuit/test_initializer.py
+++ b/test/python/circuit/test_initializer.py
@@ -25,7 +25,6 @@
     ClassicalRegister,
     BasicAer,
     transpile,
-    execute,
     assemble,
 )
 from qiskit.quantum_info import state_fidelity, Statevector, Operator
@@ -46,7 +45,8 @@ def test_uniform_superposition(self):
         qr = QuantumRegister(2, "qr")
         qc = QuantumCircuit(qr)
         qc.initialize(desired_vector, [qr[0], qr[1]])
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, desired_vector)
@@ -62,7 +62,9 @@ def test_deterministic_state(self):
         qr = QuantumRegister(2, "qr")
         qc = QuantumCircuit(qr)
         qc.initialize(desired_vector, [qr[0], qr[1]])
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, desired_vector)
@@ -87,7 +89,9 @@ def test_bell_state(self):
         qr = QuantumRegister(2, "qr")
         qc = QuantumCircuit(qr)
         qc.initialize(desired_vector, [qr[0], qr[1]])
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, desired_vector)
@@ -103,7 +107,9 @@ def test_ghz_state(self):
         qr = QuantumRegister(3, "qr")
         qc = QuantumCircuit(qr)
         qc.initialize(desired_vector, [qr[0], qr[1], qr[2]])
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, desired_vector)
@@ -120,7 +126,9 @@ def test_initialize_register(self):
         qr2 = QuantumRegister(2, "qr2")
         qc = QuantumCircuit(qr, qr2)
         qc.initialize(desired_vector, qr)
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, np.kron([1, 0, 0, 0], desired_vector))
@@ -142,7 +150,8 @@ def test_initialize_one_by_one(self):
         qc_b.initialize(qubit_0_state, [qr[0]])
         qc_b.initialize(qubit_1_state, [qr[1]])
 
-        job = execute([qc_a, qc_b], BasicAer.get_backend("statevector_simulator"))
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile([qc_a, qc_b], backend))
         result = job.result()
         statevector_a = result.get_statevector(0)
         statevector_b = result.get_statevector(1)
@@ -159,7 +168,8 @@ def test_single_qubit(self):
         qr = QuantumRegister(1, "qr")
         qc = QuantumCircuit(qr)
         qc.initialize(desired_vector, [qr[0]])
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, desired_vector)
@@ -184,7 +194,8 @@ def test_random_3qubit(self):
         qr = QuantumRegister(3, "qr")
         qc = QuantumCircuit(qr)
         qc.initialize(desired_vector, [qr[0], qr[1], qr[2]])
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, desired_vector)
@@ -217,7 +228,8 @@ def test_random_4qubit(self):
         qr = QuantumRegister(4, "qr")
         qc = QuantumCircuit(qr)
         qc.initialize(desired_vector, [qr[0], qr[1], qr[2], qr[3]])
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, desired_vector)
@@ -287,7 +299,9 @@ def test_initialize_middle_circuit(self):
         # statevector simulator does not support reset
         shots = 2000
         threshold = 0.005 * shots
-        job = execute(qc, BasicAer.get_backend("qasm_simulator"), shots=shots, seed_simulator=42)
+
+        backend = BasicAer.get_backend("qasm_simulator")
+        job = backend.run(transpile(qc, backend), shots=shots, seed_simulator=42)
         result = job.result()
         counts = result.get_counts()
         target = {"00": shots / 4, "01": shots / 4, "10": shots / 4, "11": shots / 4}
@@ -316,7 +330,8 @@ def test_math_amplitudes(self):
         qr = QuantumRegister(4, "qr")
         qc = QuantumCircuit(qr)
         qc.initialize(desired_vector, [qr[0], qr[1], qr[2], qr[3]])
-        job = execute(qc, BasicAer.get_backend("statevector_simulator"))
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc, backend))
         result = job.result()
         statevector = result.get_statevector()
         fidelity = state_fidelity(statevector, desired_vector)
@@ -338,7 +353,8 @@ def test_combiner(self):
         qc2 = QuantumCircuit(qr, cr)
         qc2.initialize(desired_vector_2, [qr[0]])
 
-        job = execute(qc1.compose(qc2), BasicAer.get_backend("statevector_simulator"))
+        backend = BasicAer.get_backend("statevector_simulator")
+        job = backend.run(transpile(qc1.compose(qc2), backend))
         result = job.result()
         quantum_state = result.get_statevector()
         fidelity = state_fidelity(quantum_state, desired_vector_2)
diff --git a/test/python/circuit/test_isometry.py b/test/python/circuit/test_isometry.py
index 3a14f7425be5..ac1b4c548ce0 100644
--- a/test/python/circuit/test_isometry.py
+++ b/test/python/circuit/test_isometry.py
@@ -21,7 +21,6 @@
 from qiskit import BasicAer
 from qiskit import QuantumCircuit
 from qiskit import QuantumRegister
-from qiskit import execute
 from qiskit.test import QiskitTestCase
 from qiskit.compiler import transpile
 from qiskit.quantum_info import Operator
@@ -66,7 +65,7 @@ def test_isometry(self, iso):
 
         # Simulate the decomposed gate
         simulator = BasicAer.get_backend("unitary_simulator")
-        result = execute(qc, simulator).result()
+        result = simulator.run(qc).result()
         unitary = result.get_unitary(qc)
         iso_from_circuit = unitary[::, 0 : 2**num_q_input]
         iso_desired = iso
@@ -108,7 +107,7 @@ def test_isometry_tolerance(self, iso):
 
         # Simulate the decomposed gate
         simulator = BasicAer.get_backend("unitary_simulator")
-        result = execute(qc, simulator).result()
+        result = simulator.run(qc).result()
         unitary = result.get_unitary(qc)
         iso_from_circuit = unitary[::, 0 : 2**num_q_input]
 
diff --git a/test/python/circuit/test_parameters.py b/test/python/circuit/test_parameters.py
index c694b39e6e7a..012498bb1a21 100644
--- a/test/python/circuit/test_parameters.py
+++ b/test/python/circuit/test_parameters.py
@@ -32,7 +32,6 @@
 from qiskit.circuit.parametertable import ParameterReferences, ParameterTable, ParameterView
 from qiskit.circuit.exceptions import CircuitError
 from qiskit.compiler import assemble, transpile
-from qiskit.execute_function import execute
 from qiskit import pulse
 from qiskit.quantum_info import Operator
 from qiskit.test import QiskitTestCase
@@ -114,7 +113,10 @@ def setUp(self):
         # TODO: delete once bind_parameters is removed from the codebase
         #  and related tests are also removed.
         warnings.filterwarnings(
-            "ignore", category=DeprecationWarning, module=r"test\.python\.circuit\.test_parameters"
+            "ignore",
+            category=DeprecationWarning,
+            module=r"test\.python\.circuit\.test_parameters",
+            message=".*bind_parameters.*",
         )
 
     def test_gate(self):
@@ -921,12 +923,8 @@ def test_transpiling_multiple_parameterized_circuits(self):
 
         circuits = [qc1, qc2]
 
-        job = execute(
-            circuits,
-            BasicAer.get_backend("unitary_simulator"),
-            shots=512,
-            parameter_binds=[{theta: 1}],
-        )
+        backend = BasicAer.get_backend("unitary_simulator")
+        job = backend.run(transpile(circuits, backend), shots=512, parameter_binds=[{theta: 1}])
 
         self.assertTrue(len(job.result().results), 2)
 
@@ -1122,7 +1120,8 @@ def test_executing_parameterized_instruction_bound_early(self, target_type):
                 bound_qc = getattr(unbound_qc, assign_fun)({theta: numpy.pi / 2})
 
                 shots = 1024
-                job = execute(bound_qc, backend=BasicAer.get_backend("qasm_simulator"), shots=shots)
+                backend = BasicAer.get_backend("qasm_simulator")
+                job = backend.run(transpile(bound_qc, backend), shots=shots)
                 self.assertDictAlmostEqual(job.result().get_counts(), {"1": shots}, 0.05 * shots)
 
     def test_num_parameters(self):
@@ -1158,7 +1157,7 @@ def test_execute_result_names(self):
 
         plist = [{theta: i} for i in range(reps)]
         simulator = BasicAer.get_backend("qasm_simulator")
-        result = execute(qc, backend=simulator, parameter_binds=plist).result()
+        result = simulator.run(transpile(qc, simulator), parameter_binds=plist).result()
         result_names = {res.name for res in result.results}
         self.assertEqual(reps, len(result_names))
 
diff --git a/test/python/circuit/test_piecewise_polynomial.py b/test/python/circuit/test_piecewise_polynomial.py
index 9bc34179ba63..639ad6f3ac6e 100644
--- a/test/python/circuit/test_piecewise_polynomial.py
+++ b/test/python/circuit/test_piecewise_polynomial.py
@@ -18,7 +18,7 @@
 from ddt import ddt, data, unpack
 
 from qiskit.test.base import QiskitTestCase
-from qiskit import BasicAer, execute
+from qiskit import BasicAer, transpile
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library.arithmetic.piecewise_polynomial_pauli_rotations import (
     PiecewisePolynomialPauliRotations,
@@ -38,7 +38,7 @@ def assertFunctionIsCorrect(self, function_circuit, reference):
         circuit.append(function_circuit.to_instruction(), list(range(circuit.num_qubits)))
 
         backend = BasicAer.get_backend("statevector_simulator")
-        statevector = execute(circuit, backend).result().get_statevector()
+        statevector = backend.run(transpile(circuit, backend)).result().get_statevector()
 
         probabilities = defaultdict(float)
         for i, statevector_amplitude in enumerate(statevector):
diff --git a/test/python/circuit/test_uc.py b/test/python/circuit/test_uc.py
index a5f0fb33e0f2..9ac5da4471dd 100644
--- a/test/python/circuit/test_uc.py
+++ b/test/python/circuit/test_uc.py
@@ -25,7 +25,7 @@
 
 from qiskit.circuit.library.generalized_gates import UCGate
 
-from qiskit import QuantumCircuit, QuantumRegister, BasicAer, execute
+from qiskit import QuantumCircuit, QuantumRegister, BasicAer
 from qiskit.test import QiskitTestCase
 from qiskit.quantum_info.random import random_unitary
 from qiskit.compiler import transpile
@@ -67,7 +67,7 @@ def test_ucg(self, squs, up_to_diagonal):
         qc = transpile(qc, basis_gates=["u1", "u3", "u2", "cx", "id"])
         # Simulate the decomposed gate
         simulator = BasicAer.get_backend("unitary_simulator")
-        result = execute(qc, simulator).result()
+        result = simulator.run(qc).result()
         unitary = result.get_unitary(qc)
         if up_to_diagonal:
             ucg = UCGate(squs, up_to_diagonal=up_to_diagonal)
@@ -86,7 +86,7 @@ def test_global_phase_ucg(self):
         qc.append(uc, q)
 
         simulator = BasicAer.get_backend("unitary_simulator")
-        result = execute(qc, simulator).result()
+        result = simulator.run(transpile(qc, simulator)).result()
         unitary = result.get_unitary(qc)
         unitary_desired = _get_ucg_matrix(gates)
 
diff --git a/test/python/circuit/test_ucx_y_z.py b/test/python/circuit/test_ucx_y_z.py
index 08f296b447c9..4450fc9e16d8 100644
--- a/test/python/circuit/test_ucx_y_z.py
+++ b/test/python/circuit/test_ucx_y_z.py
@@ -19,7 +19,7 @@
 import numpy as np
 from scipy.linalg import block_diag
 
-from qiskit import BasicAer, QuantumCircuit, QuantumRegister, execute
+from qiskit import BasicAer, QuantumCircuit, QuantumRegister
 from qiskit.test import QiskitTestCase
 
 from qiskit.quantum_info.operators.predicates import matrix_equal
@@ -67,7 +67,7 @@ def test_ucy(self, use_method):
                 qc = transpile(qc, basis_gates=["u1", "u3", "u2", "cx", "id"])
                 # Simulate the decomposed gate
                 simulator = BasicAer.get_backend("unitary_simulator")
-                result = execute(qc, simulator).result()
+                result = simulator.run(qc).result()
                 unitary = result.get_unitary(qc)
                 unitary_desired = _get_ucr_matrix(angles, rot_axis)
                 self.assertTrue(matrix_equal(unitary_desired, unitary, ignore_phase=True))
diff --git a/test/python/classical_function_compiler/test_boolean_expression.py b/test/python/classical_function_compiler/test_boolean_expression.py
index e82c1b4e8d24..12d5eebc422b 100644
--- a/test/python/classical_function_compiler/test_boolean_expression.py
+++ b/test/python/classical_function_compiler/test_boolean_expression.py
@@ -17,7 +17,7 @@
 from ddt import ddt, unpack, data
 
 from qiskit.test.base import QiskitTestCase
-from qiskit import execute, BasicAer
+from qiskit import BasicAer, transpile
 from qiskit.utils.optionals import HAS_TWEEDLEDUM
 
 if HAS_TWEEDLEDUM:
@@ -57,10 +57,11 @@ def test_synth(self, expression, expected):
         new_creg = expr_circ._create_creg(1, "c")
         expr_circ.add_register(new_creg)
         expr_circ.measure(expression.num_qubits - 1, new_creg)
+
+        backend = BasicAer.get_backend("qasm_simulator")
         [result] = (
-            execute(
-                expr_circ,
-                backend=BasicAer.get_backend("qasm_simulator"),
+            backend.run(
+                transpile(expr_circ, backend),
                 shots=1,
                 seed_simulator=14,
             )
diff --git a/test/python/compiler/test_compiler.py b/test/python/compiler/test_compiler.py
index 46f53375cd85..da697c18ea81 100644
--- a/test/python/compiler/test_compiler.py
+++ b/test/python/compiler/test_compiler.py
@@ -18,7 +18,6 @@
 from qiskit import BasicAer
 from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
 from qiskit.transpiler import PassManager
-from qiskit import execute
 from qiskit.circuit.library import U1Gate, U2Gate
 from qiskit.compiler import transpile, assemble
 from qiskit.test import QiskitTestCase
@@ -174,13 +173,16 @@ def test_example_swap_bits(self):
             qc.measure(qr0[j], ans[j])
             qc.measure(qr1[j], ans[j + n])
         # First version: no mapping
-        result = execute(
-            qc, backend=backend, coupling_map=None, shots=1024, seed_simulator=14
+        result = backend.run(
+            transpile(qc, backend), coupling_map=None, shots=1024, seed_simulator=14
         ).result()
         self.assertEqual(result.get_counts(qc), {"010000": 1024})
         # Second version: map to coupling graph
-        result = execute(
-            qc, backend=backend, coupling_map=coupling_map, shots=1024, seed_simulator=14
+        result = backend.run(
+            transpile(qc, backend, coupling_map=coupling_map),
+            coupling_map=coupling_map,
+            shots=1024,
+            seed_simulator=14,
         ).result()
         self.assertEqual(result.get_counts(qc), {"010000": 1024})
 
@@ -199,14 +201,12 @@ def test_parallel_compile(self):
         self.assertEqual(len(qobj.experiments), 10)
 
     def test_no_conflict_backend_passmanager(self):
-        """execute(qc, backend=..., passmanager=...)
-        See: https://github.com/Qiskit/qiskit-terra/issues/5037
-        """
+        """See: https://github.com/Qiskit/qiskit-terra/issues/5037"""
         backend = BasicAer.get_backend("qasm_simulator")
         qc = QuantumCircuit(2)
         qc.append(U1Gate(0), [0])
         qc.measure_all()
-        job = execute(qc, backend=backend, pass_manager=PassManager())
+        job = backend.run(PassManager().run(qc))
         result = job.result().get_counts()
         self.assertEqual(result, {"00": 1024})
 
@@ -302,21 +302,20 @@ def test_mapper_overoptimization(self):
         coupling_map = [[0, 2], [1, 2], [2, 3]]
         shots = 1000
 
-        result1 = execute(
-            circ,
-            backend=self.backend,
-            coupling_map=coupling_map,
+        result1 = self.backend.run(
+            transpile(circ, backend=self.backend, coupling_map=coupling_map, seed_transpiler=8),
             seed_simulator=self.seed_simulator,
-            seed_transpiler=8,
             shots=shots,
         )
         count1 = result1.result().get_counts()
-        result2 = execute(
-            circ,
-            backend=self.backend,
-            coupling_map=None,
+        result2 = self.backend.run(
+            transpile(
+                circ,
+                backend=self.backend,
+                coupling_map=None,
+                seed_transpiler=8,
+            ),
             seed_simulator=self.seed_simulator,
-            seed_transpiler=8,
             shots=shots,
         )
         count2 = result2.result().get_counts()
@@ -376,9 +375,12 @@ def test_grovers_circuit(self):
         circuit.measure(qr[0], cr[0])
         circuit.measure(qr[1], cr[1])
 
-        result = execute(
-            circuit,
-            backend=self.backend,
+        result = self.backend.run(
+            transpile(
+                circuit,
+                self.backend,
+                coupling_map=coupling_map,
+            ),
             coupling_map=coupling_map,
             seed_simulator=self.seed_simulator,
             shots=shots,
@@ -427,9 +429,11 @@ def test_math_domain_error(self):
 
         coupling_map = [[0, 2], [1, 2], [2, 3]]
         shots = 2000
-        job = execute(
-            circ,
-            backend=self.backend,
+        job = self.backend.run(
+            transpile(
+                circ,
+                backend=self.backend,
+            ),
             coupling_map=coupling_map,
             seed_simulator=self.seed_simulator,
             shots=shots,
@@ -445,10 +449,8 @@ def test_random_parameter_circuit(self):
         circ = QuantumCircuit.from_qasm_file(os.path.join(qasm_dir, "random_n5_d5.qasm"))
         coupling_map = [[0, 1], [1, 2], [2, 3], [3, 4]]
         shots = 1024
-        qobj = execute(
-            circ,
-            backend=self.backend,
-            coupling_map=coupling_map,
+        qobj = self.backend.run(
+            transpile(circ, backend=self.backend, coupling_map=coupling_map, seed_transpiler=42),
             shots=shots,
             seed_simulator=self.seed_simulator,
         )
diff --git a/test/python/opflow/test_state_construction.py b/test/python/opflow/test_state_construction.py
index 05d595c18da4..97c24e94adda 100644
--- a/test/python/opflow/test_state_construction.py
+++ b/test/python/opflow/test_state_construction.py
@@ -16,7 +16,7 @@
 from test.python.opflow import QiskitOpflowTestCase
 import numpy as np
 
-from qiskit import QuantumCircuit, BasicAer, execute
+from qiskit import QuantumCircuit, BasicAer
 from qiskit.circuit import ParameterVector
 from qiskit.quantum_info import Statevector
 
@@ -85,13 +85,11 @@ def test_qiskit_result_instantiation(self):
         qc = QuantumCircuit(3)
         # REMEMBER: This is Qubit 2 in Operator land.
         qc.h(0)
-        sv_res = execute(qc, BasicAer.get_backend("statevector_simulator")).result()
+        sv_res = BasicAer.get_backend("statevector_simulator").run(qc).result()
         sv_vector = sv_res.get_statevector()
         qc_op = PrimitiveOp(qc) @ Zero
 
-        qasm_res = execute(
-            qc_op.to_circuit(meas=True), BasicAer.get_backend("qasm_simulator")
-        ).result()
+        qasm_res = BasicAer.get_backend("qasm_simulator").run(qc_op.to_circuit(meas=True)).result()
 
         np.testing.assert_array_almost_equal(
             StateFn(sv_res).to_matrix(), [0.5**0.5, 0.5**0.5, 0, 0, 0, 0, 0, 0]
diff --git a/test/python/providers/test_fake_backends.py b/test/python/providers/test_fake_backends.py
index 078a5110f818..eaaf79bfc462 100644
--- a/test/python/providers/test_fake_backends.py
+++ b/test/python/providers/test_fake_backends.py
@@ -25,7 +25,6 @@
 from qiskit.compiler import assemble
 from qiskit.compiler import transpile
 from qiskit.exceptions import QiskitError
-from qiskit.execute_function import execute
 from qiskit.test.base import QiskitTestCase
 from qiskit.providers.fake_provider import (
     FakeProviderForBackendV2,
@@ -102,12 +101,10 @@ def setUpClass(cls):
     def test_circuit_on_fake_backend_v2(self, backend, optimization_level):
         if not optionals.HAS_AER and backend.num_qubits > 20:
             self.skipTest("Unable to run fake_backend %s without qiskit-aer" % backend.backend_name)
-        job = execute(
-            self.circuit,
-            backend,
+        job = backend.run(
+            transpile(self.circuit, backend, seed_transpiler=42),
             optimization_level=optimization_level,
             seed_simulator=42,
-            seed_transpiler=42,
         )
         result = job.result()
         counts = result.get_counts()
@@ -126,12 +123,10 @@ def test_circuit_on_fake_backend(self, backend, optimization_level):
                 "Unable to run fake_backend %s without qiskit-aer"
                 % backend.configuration().backend_name
             )
-        job = execute(
-            self.circuit,
-            backend,
+        job = backend.run(
+            transpile(self.circuit, backend, seed_transpiler=42),
             optimization_level=optimization_level,
             seed_simulator=42,
-            seed_transpiler=42,
         )
         result = job.result()
         counts = result.get_counts()
diff --git a/test/python/quantum_info/test_analyzation.py b/test/python/quantum_info/test_analyzation.py
index 9feed4f390f1..fe44862402cb 100644
--- a/test/python/quantum_info/test_analyzation.py
+++ b/test/python/quantum_info/test_analyzation.py
@@ -14,8 +14,7 @@
 
 import unittest
 
-import qiskit
-from qiskit import BasicAer
+from qiskit import BasicAer, QuantumCircuit, ClassicalRegister, QuantumRegister, transpile
 from qiskit.quantum_info.analysis.average import average_data
 from qiskit.quantum_info.analysis.make_observable import make_dict_observable
 from qiskit.quantum_info.analysis import hellinger_fidelity
@@ -27,16 +26,16 @@ class TestAnalyzation(QiskitTestCase):
 
     def test_average_data_dict_observable(self):
         """Test average_data for dictionary observable input"""
-        qr = qiskit.QuantumRegister(2)
-        cr = qiskit.ClassicalRegister(2)
-        qc = qiskit.QuantumCircuit(qr, cr, name="qc")
+        qr = QuantumRegister(2)
+        cr = ClassicalRegister(2)
+        qc = QuantumCircuit(qr, cr, name="qc")
         qc.h(qr[0])
         qc.cx(qr[0], qr[1])
         qc.measure(qr[0], cr[0])
         qc.measure(qr[1], cr[1])
         shots = 10000
         backend = BasicAer.get_backend("qasm_simulator")
-        result = qiskit.execute(qc, backend, shots=shots).result()
+        result = backend.run(qc, shots=shots).result()
         counts = result.get_counts(qc)
         observable = {"00": 1, "11": 1, "01": -1, "10": -1}
         mean_zz = average_data(counts=counts, observable=observable)
@@ -50,9 +49,9 @@ def test_average_data_dict_observable(self):
 
     def test_average_data_list_observable(self):
         """Test average_data for list observable input."""
-        qr = qiskit.QuantumRegister(3)
-        cr = qiskit.ClassicalRegister(3)
-        qc = qiskit.QuantumCircuit(qr, cr, name="qc")
+        qr = QuantumRegister(3)
+        cr = ClassicalRegister(3)
+        qc = QuantumCircuit(qr, cr, name="qc")
         qc.h(qr[0])
         qc.cx(qr[0], qr[1])
         qc.cx(qr[0], qr[2])
@@ -61,7 +60,7 @@ def test_average_data_list_observable(self):
         qc.measure(qr[2], cr[2])
         shots = 10000
         backend = BasicAer.get_backend("qasm_simulator")
-        result = qiskit.execute(qc, backend, shots=shots).result()
+        result = backend.run(qc, shots=shots).result()
         counts = result.get_counts(qc)
         observable = [1, -1, -1, 1, -1, 1, 1, -1]
         mean_zzz = average_data(counts=counts, observable=observable)
@@ -78,16 +77,16 @@ def test_average_data_list_observable(self):
 
     def test_average_data_matrix_observable(self):
         """Test average_data for matrix observable input."""
-        qr = qiskit.QuantumRegister(2)
-        cr = qiskit.ClassicalRegister(2)
-        qc = qiskit.QuantumCircuit(qr, cr, name="qc")
+        qr = QuantumRegister(2)
+        cr = ClassicalRegister(2)
+        qc = QuantumCircuit(qr, cr, name="qc")
         qc.h(qr[0])
         qc.cx(qr[0], qr[1])
         qc.measure(qr[0], cr[0])
         qc.measure(qr[1], cr[1])
         shots = 10000
         backend = BasicAer.get_backend("qasm_simulator")
-        result = qiskit.execute(qc, backend, shots=shots).result()
+        result = backend.run(qc, shots=shots).result()
         counts = result.get_counts(qc)
         observable = [[1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]]
         mean_zz = average_data(counts=counts, observable=observable)
@@ -125,7 +124,7 @@ def test_make_dict_observable(self):
 
     def test_hellinger_fidelity_same(self):
         """Test hellinger fidelity is one for same dist."""
-        qc = qiskit.QuantumCircuit(5, 5)
+        qc = QuantumCircuit(5, 5)
         qc.h(2)
         qc.cx(2, 1)
         qc.cx(2, 3)
@@ -135,7 +134,7 @@ def test_hellinger_fidelity_same(self):
 
         sim = BasicAer.get_backend("qasm_simulator")
 
-        res = qiskit.execute(qc, sim).result()
+        res = sim.run(qc).result()
 
         ans = hellinger_fidelity(res.get_counts(), res.get_counts())
 
@@ -157,7 +156,7 @@ def test_hellinger_fidelity_no_overlap(self):
         #                     └───┘ ║  ║  ║  ║ └╥┘
         # c: 5/═════════════════════╩══╩══╩══╩══╩═
         #                           0  1  2  3  4
-        qc = qiskit.QuantumCircuit(5, 5)
+        qc = QuantumCircuit(5, 5)
         qc.h(2)
         qc.cx(2, 1)
         qc.cx(2, 3)
@@ -178,7 +177,7 @@ def test_hellinger_fidelity_no_overlap(self):
         #       └╥┘                 ║  ║  ║  ║
         # c: 5/══╩══════════════════╩══╩══╩══╩═
         #        4                  0  1  2  3
-        qc2 = qiskit.QuantumCircuit(5, 5)
+        qc2 = QuantumCircuit(5, 5)
         qc2.h(2)
         qc2.cx(2, 1)
         qc2.y(2)
@@ -188,8 +187,8 @@ def test_hellinger_fidelity_no_overlap(self):
 
         sim = BasicAer.get_backend("qasm_simulator")
 
-        res1 = qiskit.execute(qc, sim).result()
-        res2 = qiskit.execute(qc2, sim).result()
+        res1 = sim.run(qc).result()
+        res2 = sim.run(transpile(qc2, sim)).result()
 
         ans = hellinger_fidelity(res1.get_counts(), res2.get_counts())
 
diff --git a/test/python/quantum_info/test_local_invariance.py b/test/python/quantum_info/test_local_invariance.py
index 3202e799b646..320a24b06221 100644
--- a/test/python/quantum_info/test_local_invariance.py
+++ b/test/python/quantum_info/test_local_invariance.py
@@ -16,8 +16,7 @@
 import unittest
 
 from numpy.testing import assert_allclose
-from qiskit.execute_function import execute
-from qiskit.circuit import QuantumCircuit, QuantumRegister
+from qiskit import QuantumCircuit, QuantumRegister, transpile
 from qiskit.test import QiskitTestCase
 from qiskit.providers.basicaer import UnitarySimulatorPy
 from qiskit.quantum_info.synthesis.local_invariance import two_qubit_local_invariants
@@ -34,14 +33,14 @@ def test_2q_local_invariance_simple(self):
 
         qr = QuantumRegister(2, name="q")
         qc = QuantumCircuit(qr)
-        U = execute(qc, sim).result().get_unitary()
+        U = sim.run(qc).result().get_unitary()
         vec = two_qubit_local_invariants(U)
         assert_allclose(vec, [1, 0, 3])
 
         qr = QuantumRegister(2, name="q")
         qc = QuantumCircuit(qr)
         qc.cx(qr[1], qr[0])
-        U = execute(qc, sim).result().get_unitary()
+        U = sim.run(qc).result().get_unitary()
         vec = two_qubit_local_invariants(U)
         assert_allclose(vec, [0, 0, 1])
 
@@ -49,14 +48,14 @@ def test_2q_local_invariance_simple(self):
         qc = QuantumCircuit(qr)
         qc.cx(qr[1], qr[0])
         qc.cx(qr[0], qr[1])
-        U = execute(qc, sim).result().get_unitary()
+        U = sim.run(qc).result().get_unitary()
         vec = two_qubit_local_invariants(U)
         assert_allclose(vec, [0, 0, -1])
 
         qr = QuantumRegister(2, name="q")
         qc = QuantumCircuit(qr)
         qc.swap(qr[1], qr[0])
-        U = execute(qc, sim).result().get_unitary()
+        U = sim.run(transpile(qc, sim)).result().get_unitary()
         vec = two_qubit_local_invariants(U)
         assert_allclose(vec, [-1, 0, -3])
 
diff --git a/test/python/quantum_info/test_synthesis.py b/test/python/quantum_info/test_synthesis.py
index 2d134967e34f..4bbf8304a126 100644
--- a/test/python/quantum_info/test_synthesis.py
+++ b/test/python/quantum_info/test_synthesis.py
@@ -24,7 +24,7 @@
 import scipy.stats
 from ddt import ddt, data
 
-from qiskit import execute, QiskitError, transpile
+from qiskit import QiskitError, transpile
 from qiskit.circuit import QuantumCircuit, QuantumRegister
 from qiskit.converters import dag_to_circuit, circuit_to_dag
 from qiskit.circuit.library import (
@@ -52,7 +52,6 @@
     RZGate,
     UnitaryGate,
 )
-from qiskit.providers.basicaer import UnitarySimulatorPy
 from qiskit.quantum_info.operators import Operator
 from qiskit.quantum_info.random import random_unitary
 from qiskit.quantum_info.synthesis.one_qubit_decompose import OneQubitEulerDecomposer
@@ -270,8 +269,7 @@ def check_exact_decomposition(
         decomp_circuit = decomposer(target_unitary, _num_basis_uses=num_basis_uses)
         if num_basis_uses is not None:
             self.assertEqual(num_basis_uses, decomp_circuit.count_ops().get("unitary", 0))
-        result = execute(decomp_circuit, UnitarySimulatorPy(), optimization_level=0).result()
-        decomp_unitary = result.get_unitary()
+        decomp_unitary = Operator(decomp_circuit).data
         maxdist = np.max(np.abs(target_unitary - decomp_unitary))
         self.assertTrue(
             np.abs(maxdist) < tolerance,
@@ -1100,10 +1098,7 @@ def test_approx_supercontrolled_decompose_random(self, seed):
         for i in range(4):
             with self.subTest(i=i):
                 decomp_circuit = decomposer(tgt, _num_basis_uses=i)
-                result = execute(
-                    decomp_circuit, UnitarySimulatorPy(), optimization_level=0
-                ).result()
-                decomp_unitary = result.get_unitary()
+                decomp_unitary = Operator(decomp_circuit).data
                 tr_actual = np.trace(decomp_unitary.conj().T @ tgt)
                 self.assertAlmostEqual(
                     traces_pred[i],
@@ -1123,8 +1118,7 @@ def test_cx_equivalence_0cx(self, seed=0):
         qc.u(rnd[0], rnd[1], rnd[2], qr[0])
         qc.u(rnd[3], rnd[4], rnd[5], qr[1])
 
-        sim = UnitarySimulatorPy()
-        unitary = execute(qc, sim).result().get_unitary()
+        unitary = Operator(qc).data
         self.assertEqual(two_qubit_cnot_decompose.num_basis_gates(unitary), 0)
         self.assertTrue(Operator(two_qubit_cnot_decompose(unitary)).equiv(unitary))
 
@@ -1144,8 +1138,7 @@ def test_cx_equivalence_1cx(self, seed=1):
         qc.u(rnd[6], rnd[7], rnd[8], qr[0])
         qc.u(rnd[9], rnd[10], rnd[11], qr[1])
 
-        sim = UnitarySimulatorPy()
-        unitary = execute(qc, sim).result().get_unitary()
+        unitary = Operator(qc).data
         self.assertEqual(two_qubit_cnot_decompose.num_basis_gates(unitary), 1)
         self.assertTrue(Operator(two_qubit_cnot_decompose(unitary)).equiv(unitary))
 
@@ -1170,8 +1163,7 @@ def test_cx_equivalence_2cx(self, seed=2):
         qc.u(rnd[12], rnd[13], rnd[14], qr[0])
         qc.u(rnd[15], rnd[16], rnd[17], qr[1])
 
-        sim = UnitarySimulatorPy()
-        unitary = execute(qc, sim).result().get_unitary()
+        unitary = Operator(qc).data
         self.assertEqual(two_qubit_cnot_decompose.num_basis_gates(unitary), 2)
         self.assertTrue(Operator(two_qubit_cnot_decompose(unitary)).equiv(unitary))
 
@@ -1201,8 +1193,7 @@ def test_cx_equivalence_3cx(self, seed=3):
         qc.u(rnd[18], rnd[19], rnd[20], qr[0])
         qc.u(rnd[21], rnd[22], rnd[23], qr[1])
 
-        sim = UnitarySimulatorPy()
-        unitary = execute(qc, sim).result().get_unitary()
+        unitary = Operator(qc).data
         self.assertEqual(two_qubit_cnot_decompose.num_basis_gates(unitary), 3)
         self.assertTrue(Operator(two_qubit_cnot_decompose(unitary)).equiv(unitary))
 
diff --git a/test/python/tools/monitor/test_job_monitor.py b/test/python/tools/monitor/test_job_monitor.py
index 17409ac1022a..95a93e65ea76 100644
--- a/test/python/tools/monitor/test_job_monitor.py
+++ b/test/python/tools/monitor/test_job_monitor.py
@@ -33,7 +33,8 @@ def test_job_monitor(self):
         qc.cx(qreg[0], qreg[1])
         qc.measure(qreg, creg)
         backend = BasicAer.get_backend("qasm_simulator")
-        job_sim = execute([qc] * 10, backend)
+        with self.assertWarns(DeprecationWarning):
+            job_sim = execute([qc] * 10, backend)
         output = io.StringIO()
         job_monitor(job_sim, output=output)
         self.assertEqual(job_sim.status().name, "DONE")
diff --git a/test/python/utils/mitigation/test_meas.py b/test/python/utils/mitigation/test_meas.py
index a5105b5ba4d5..d91be70c5de5 100644
--- a/test/python/utils/mitigation/test_meas.py
+++ b/test/python/utils/mitigation/test_meas.py
@@ -286,7 +286,8 @@ def test_ideal_meas_cal(self):
 
                 # Perform an ideal execution on the generated circuits
                 backend = AerSimulator()
-                job = qiskit.execute(meas_calibs, backend=backend, shots=self.shots)
+                with self.assertWarns(DeprecationWarning):
+                    job = qiskit.execute(meas_calibs, backend=backend, shots=self.shots)
                 cal_results = job.result()
 
                 with self.assertWarns(DeprecationWarning):
@@ -339,13 +340,14 @@ def test_meas_cal_on_circuit(self):
 
         # Run the calibration circuits
         backend = AerSimulator()
-        job = qiskit.execute(
-            meas_calibs,
-            backend=backend,
-            shots=self.shots,
-            seed_simulator=SEED,
-            seed_transpiler=SEED,
-        )
+        with self.assertWarns(DeprecationWarning):
+            job = qiskit.execute(
+                meas_calibs,
+                backend=backend,
+                shots=self.shots,
+                seed_simulator=SEED,
+                seed_transpiler=SEED,
+            )
         cal_results = job.result()
 
         with self.assertWarns(DeprecationWarning):
@@ -353,10 +355,10 @@ def test_meas_cal_on_circuit(self):
             meas_cal = CompleteMeasFitter(cal_results, state_labels)
         # Calculate the fidelity
         fidelity = meas_cal.readout_fidelity()
-
-        job = qiskit.execute(
-            [ghz], backend=backend, shots=self.shots, seed_simulator=SEED, seed_transpiler=SEED
-        )
+        with self.assertWarns(DeprecationWarning):
+            job = qiskit.execute(
+                [ghz], backend=backend, shots=self.shots, seed_simulator=SEED, seed_transpiler=SEED
+            )
         results = job.result()
 
         # Predicted equally distributed results
@@ -403,7 +405,9 @@ def test_ideal_tensored_meas_cal(self):
 
         # Perform an ideal execution on the generated circuits
         backend = AerSimulator()
-        cal_results = qiskit.execute(meas_calibs, backend=backend, shots=self.shots).result()
+        with self.assertWarns(DeprecationWarning):
+            job = qiskit.execute(meas_calibs, backend=backend, shots=self.shots)
+        cal_results = job.result()
 
         with self.assertWarns(DeprecationWarning):
             # Make calibration matrices
@@ -459,13 +463,15 @@ def test_tensored_meas_cal_on_circuit(self):
 
         # Run the calibration circuits
         backend = AerSimulator()
-        cal_results = qiskit.execute(
-            meas_calibs,
-            backend=backend,
-            shots=self.shots,
-            seed_simulator=SEED,
-            seed_transpiler=SEED,
-        ).result()
+        with self.assertWarns(DeprecationWarning):
+            job = qiskit.execute(
+                meas_calibs,
+                backend=backend,
+                shots=self.shots,
+                seed_simulator=SEED,
+                seed_transpiler=SEED,
+            )
+        cal_results = job.result()
 
         with self.assertWarns(DeprecationWarning):
             # Make a calibration matrix
@@ -473,9 +479,10 @@ def test_tensored_meas_cal_on_circuit(self):
         # Calculate the fidelity
         fidelity = meas_cal.readout_fidelity(0) * meas_cal.readout_fidelity(1)
 
-        results = qiskit.execute(
-            [ghz], backend=backend, shots=self.shots, seed_simulator=SEED, seed_transpiler=SEED
-        ).result()
+        with self.assertWarns(DeprecationWarning):
+            results = qiskit.execute(
+                [ghz], backend=backend, shots=self.shots, seed_simulator=SEED, seed_transpiler=SEED
+            ).result()
 
         # Predicted equally distributed results
         predicted_results = {"000": 0.5, "111": 0.5}
diff --git a/test/randomized/test_synthesis.py b/test/randomized/test_synthesis.py
index 2b3e387b9baf..ef7ac3304f09 100644
--- a/test/randomized/test_synthesis.py
+++ b/test/randomized/test_synthesis.py
@@ -16,7 +16,6 @@
 from hypothesis import given, strategies, settings
 import numpy as np
 
-from qiskit import execute
 from qiskit.circuit import QuantumCircuit, QuantumRegister
 from qiskit.circuit.library import UnitaryGate
 from qiskit.providers.basicaer import UnitarySimulatorPy
@@ -74,7 +73,7 @@ def test_cx_equivalence_0cx_random(self, rnd, seed):
         qc.u(rnd[3], rnd[4], rnd[5], qr[1])
 
         sim = UnitarySimulatorPy()
-        unitary = execute(qc, sim, seed_simulator=seed).result().get_unitary()
+        unitary = sim.run(qc, seed_simulator=seed).result().get_unitary()
         self.assertEqual(two_qubit_cnot_decompose.num_basis_gates(unitary), 0)
 
     @given(strategies.tuples(*[rotation] * 12), seed)
@@ -92,7 +91,7 @@ def test_cx_equivalence_1cx_random(self, rnd, seed):
         qc.u(rnd[9], rnd[10], rnd[11], qr[1])
 
         sim = UnitarySimulatorPy()
-        unitary = execute(qc, sim, seed_simulator=seed).result().get_unitary()
+        unitary = sim.run(qc, seed_simulator=seed).result().get_unitary()
         self.assertEqual(two_qubit_cnot_decompose.num_basis_gates(unitary), 1)
 
     @given(strategies.tuples(*[rotation] * 18), seed)
@@ -115,7 +114,7 @@ def test_cx_equivalence_2cx_random(self, rnd, seed):
         qc.u(rnd[15], rnd[16], rnd[17], qr[1])
 
         sim = UnitarySimulatorPy()
-        unitary = execute(qc, sim, seed_simulator=seed).result().get_unitary()
+        unitary = sim.run(qc, seed_simulator=seed).result().get_unitary()
         self.assertEqual(two_qubit_cnot_decompose.num_basis_gates(unitary), 2)
 
     @given(strategies.tuples(*[rotation] * 24), seed)
@@ -143,7 +142,7 @@ def test_cx_equivalence_3cx_random(self, rnd, seed):
         qc.u(rnd[21], rnd[22], rnd[23], qr[1])
 
         sim = UnitarySimulatorPy()
-        unitary = execute(qc, sim, seed_simulator=seed).result().get_unitary()
+        unitary = sim.run(qc, seed_simulator=seed).result().get_unitary()
         self.assertEqual(two_qubit_cnot_decompose.num_basis_gates(unitary), 3)
 
 
diff --git a/test/visual/mpl/graph/test_graph_matplotlib_drawer.py b/test/visual/mpl/graph/test_graph_matplotlib_drawer.py
index 46e9c4b39653..fe3a6260d83f 100644
--- a/test/visual/mpl/graph/test_graph_matplotlib_drawer.py
+++ b/test/visual/mpl/graph/test_graph_matplotlib_drawer.py
@@ -18,7 +18,7 @@
 from contextlib import contextmanager
 from pathlib import Path
 
-from qiskit import BasicAer, execute
+from qiskit import BasicAer, transpile
 from qiskit.test import QiskitTestCase
 from qiskit import QuantumCircuit
 from qiskit.utils import optionals
@@ -102,7 +102,7 @@ def test_plot_bloch_multivector(self):
 
         # getting the state using backend
         backend = BasicAer.get_backend("statevector_simulator")
-        result = execute(circuit, backend).result()
+        result = backend.run(circuit).result()
         state = result.get_statevector(circuit)
 
         fname = "bloch_multivector.png"
@@ -124,7 +124,7 @@ def test_plot_state_hinton(self):
 
         # getting the state using backend
         backend = BasicAer.get_backend("statevector_simulator")
-        result = execute(circuit, backend).result()
+        result = backend.run(circuit).result()
         state = result.get_statevector(circuit)
 
         fname = "hinton.png"
@@ -146,7 +146,7 @@ def test_plot_state_qsphere(self):
 
         # getting the state using backend
         backend = BasicAer.get_backend("statevector_simulator")
-        result = execute(circuit, backend).result()
+        result = backend.run(circuit).result()
         state = result.get_statevector(circuit)
 
         fname = "qsphere.png"
@@ -168,7 +168,7 @@ def test_plot_state_city(self):
 
         # getting the state using backend
         backend = BasicAer.get_backend("statevector_simulator")
-        result = execute(circuit, backend).result()
+        result = backend.run(circuit).result()
         state = result.get_statevector(circuit)
 
         fname = "state_city.png"
@@ -190,7 +190,7 @@ def test_plot_state_paulivec(self):
 
         # getting the state using backend
         backend = BasicAer.get_backend("statevector_simulator")
-        result = execute(circuit, backend).result()
+        result = backend.run(circuit).result()
         state = result.get_statevector(circuit)
 
         fname = "paulivec.png"
@@ -648,7 +648,7 @@ def test_plot_bloch_multivector_figsize_improvements(self):
 
         # getting the state using backend
         backend = BasicAer.get_backend("statevector_simulator")
-        result = execute(circuit, backend).result()
+        result = backend.run(transpile(circuit, backend)).result()
         state = result.get_statevector(circuit)
 
         fname = "bloch_multivector_figsize_improvements.png"