diff --git a/pennylane_qiskit/qiskit_device.py b/pennylane_qiskit/qiskit_device.py
index 09b0a7588..bdf959695 100644
--- a/pennylane_qiskit/qiskit_device.py
+++ b/pennylane_qiskit/qiskit_device.py
@@ -1,4 +1,4 @@
-# Copyright 2019 Xanadu Quantum Technologies Inc.
+# Copyright 2019-2021 Xanadu Quantum Technologies Inc.
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -203,7 +203,17 @@ def reset(self):
         self._current_job = None
         self._state = None  # statevector of a simulator backend
 
-    def apply(self, operations, **kwargs):
+    def create_circuit_object(self, operations, **kwargs):
+        """Builds the circuit objects based on the operations and measurements
+        specified to apply.
+
+        Args:
+            operations (list[~.Operation]): operations to apply to the device
+
+        Keyword args:
+            rotations (list[~.Operation]): Operations that rotate the circuit
+                pre-measurement into the eigenbasis of the observables.
+        """
         rotations = kwargs.get("rotations", [])
 
         applied_operations = self.apply_operations(operations)
@@ -222,6 +232,10 @@ def apply(self, operations, **kwargs):
         elif "aer" in self.backend_name:
             self._circuit.save_state()
 
+    def apply(self, operations, **kwargs):
+
+        self.create_circuit_object(operations, **kwargs)
+
         # These operations need to run for all devices
         compiled_circuit = self.compile()
         self.run(compiled_circuit)
@@ -305,20 +319,21 @@ def run(self, qcirc):
         if self.backend_name in self._state_backends:
             self._state = self._get_state(result)
 
-    def _get_state(self, result):
+    def _get_state(self, result, experiment=None):
         """Returns the statevector for state simulator backends.
 
         Args:
             result (qiskit.Result): result object
+            experiment (str): the name of the experiment to get the state for
 
         Returns:
             array[float]: size ``(2**num_wires,)`` statevector
         """
         if "statevector" in self.backend_name:
-            state = np.asarray(result.get_statevector())
+            state = np.asarray(result.get_statevector(experiment))
 
         elif "unitary" in self.backend_name:
-            unitary = np.asarray(result.get_unitary())
+            unitary = np.asarray(result.get_unitary(experiment))
             initial_state = np.zeros([2 ** self.num_wires])
             initial_state[0] = 1
 
@@ -327,7 +342,18 @@ def _get_state(self, result):
         # reverse qubit order to match PennyLane convention
         return state.reshape([2] * self.num_wires).T.flatten()
 
-    def generate_samples(self):
+    def generate_samples(self, circuit=None):
+        r"""Returns the computational basis samples generated for all wires.
+
+        Note that PennyLane uses the convention :math:`|q_0,q_1,\dots,q_{N-1}\rangle` where
+        :math:`q_0` is the most significant bit.
+
+        Args:
+            circuit (str): the name of the circuit to get the state for
+
+        Returns:
+             array[complex]: array of samples in the shape ``(dev.shots, dev.num_wires)``
+        """
 
         # branch out depending on the type of backend
         if self.backend_name in self._state_backends:
@@ -335,7 +361,7 @@ def generate_samples(self):
             return super().generate_samples()
 
         # hardware or hardware simulator
-        samples = self._current_job.result().get_memory()
+        samples = self._current_job.result().get_memory(circuit)
 
         # reverse qubit order to match PennyLane convention
         return np.vstack([np.array([int(i) for i in s[::-1]]) for s in samples])
@@ -350,3 +376,39 @@ def analytic_probability(self, wires=None):
 
         prob = self.marginal_prob(np.abs(self._state) ** 2, wires)
         return prob
+
+    def batch_execute(self, circuits):
+
+        compiled_circuits = []
+
+        # Compile each circuit object
+        for circuit in circuits:
+
+            # We need to reset the device here, else it will
+            # not start the next computation in the zero state
+            self.reset()
+            self.create_circuit_object(circuit.operations, rotations=circuit.diagonalizing_gates)
+
+            compiled_circ = self.compile()
+            compiled_circ.name = f"circ{len(compiled_circuits)}"
+            compiled_circuits.append(compiled_circ)
+
+        # Send the batch of circuit objects using backend.run
+        self._current_job = self.backend.run(compiled_circuits, shots=self.shots, **self.run_args)
+        result = self._current_job.result()
+
+        # Compute statistics using the state and/or samples
+        results = []
+        for circuit, circuit_obj in zip(circuits, compiled_circuits):
+
+            if self.backend_name in self._state_backends:
+                self._state = self._get_state(result, experiment=circuit_obj)
+
+            # generate computational basis samples
+            if self.shots is not None or circuit.is_sampled:
+                self._samples = self.generate_samples(circuit_obj)
+
+            res = self.statistics(circuit.observables)
+            results.append(res)
+
+        return results
diff --git a/tests/test_ibmq.py b/tests/test_ibmq.py
index c3ab7558d..d9e1c83d5 100644
--- a/tests/test_ibmq.py
+++ b/tests/test_ibmq.py
@@ -14,6 +14,8 @@
 
 @pytest.fixture
 def token():
+    """A fixture loading the IBMQ token from the IBMQX_TOKEN_TEST environment
+    variable."""
     t = os.getenv("IBMQX_TOKEN_TEST", None)
 
     if t is None:
@@ -131,6 +133,7 @@ def test_custom_provider_hub_group_project(monkeypatch):
 
 
 def test_load_from_disk(token):
+    """Test loading the account credentials and the device from disk."""
     IBMQ.save_account(token)
     dev = IBMQDevice(wires=1)
     assert dev.provider.credentials.is_ibmq()
@@ -138,6 +141,7 @@ def test_load_from_disk(token):
 
 
 def test_account_error(monkeypatch):
+    """Test that an error is raised if there is no active IBMQ account."""
 
     # Token is passed such that the test is skipped if no token was provided
     with pytest.raises(IBMQAccountError, match="No active IBM Q account"):
@@ -148,13 +152,14 @@ def test_account_error(monkeypatch):
 
 @pytest.mark.parametrize("shots", [1000])
 def test_simple_circuit(token, tol, shots):
+    """Test executing a simple circuit submitted to IBMQ."""
     IBMQ.enable_account(token)
     dev = IBMQDevice(wires=2, backend="ibmq_qasm_simulator", shots=shots)
 
     @qml.qnode(dev)
     def circuit(theta, phi):
         qml.RX(theta, wires=0)
-        qml.RX(phi, wires=0)
+        qml.RX(phi, wires=1)
         qml.CNOT(wires=[0, 1])
         return qml.expval(qml.PauliZ(0)), qml.expval(qml.PauliZ(1))
 
@@ -166,6 +171,70 @@ def circuit(theta, phi):
     assert np.allclose(res, expected, **tol)
 
 
+@pytest.mark.parametrize("shots", [1000])
+def test_simple_circuit_with_batch_params(token, tol, shots, mocker):
+    """Test that executing a simple circuit with batched parameters is
+    submitted to IBMQ once."""
+    IBMQ.enable_account(token)
+    dev = IBMQDevice(wires=2, backend="ibmq_qasm_simulator", shots=shots)
+
+    @qml.batch_params
+    @qml.qnode(dev)
+    def circuit(theta, phi):
+        qml.RX(theta, wires=0)
+        qml.RX(phi, wires=1)
+        qml.CNOT(wires=[0, 1])
+        return qml.expval(qml.PauliZ(0)), qml.expval(qml.PauliZ(1))
+
+    # Check that we run only once
+    spy1 = mocker.spy(dev, "batch_execute")
+    spy2 = mocker.spy(dev.backend, "run")
+
+    # Batch the input parameters
+    batch_dim = 3
+    theta = np.linspace(0, 0.543, batch_dim)
+    phi = np.linspace(0, 0.123, batch_dim)
+
+    res = circuit(theta, phi)
+    assert np.allclose(res[:, 0], np.cos(theta), **tol)
+    assert np.allclose(res[:, 1], np.cos(theta) * np.cos(phi), **tol)
+
+    # Check that IBMQBackend.run was called once
+    assert spy1.call_count == 1
+    assert spy2.call_count == 1
+
+
+@pytest.mark.parametrize("shots", [1000])
+def test_batch_execute_parameter_shift(token, tol, shots, mocker):
+    """Test that devices provide correct result computing the gradient of a
+    circuit using the parameter-shift rule and the batch execution pipeline."""
+    IBMQ.enable_account(token)
+    dev = IBMQDevice(wires=3, backend="ibmq_qasm_simulator", shots=shots)
+
+    spy1 = mocker.spy(dev, "batch_execute")
+    spy2 = mocker.spy(dev.backend, "run")
+
+    @qml.qnode(dev, diff_method="parameter-shift")
+    def circuit(x, y):
+        qml.RX(x, wires=[0])
+        qml.RY(y, wires=[1])
+        qml.CNOT(wires=[0, 1])
+        return qml.expval(qml.PauliZ(0) @ qml.PauliX(1) @ qml.PauliZ(2))
+
+    x = qml.numpy.array(0.543, requires_grad=True)
+    y = qml.numpy.array(0.123, requires_grad=True)
+
+    res = qml.grad(circuit)(x,y)
+    expected = np.array([[-np.sin(y) * np.sin(x), np.cos(y) * np.cos(x)]])
+    assert np.allclose(res, expected, **tol)
+
+    # Check that QiskitDevice.batch_execute was called once
+    assert spy1.call_count == 1
+
+    # Check that run was called twice: for the partial derivatives and for
+    # running the circuit
+    assert spy2.call_count == 2
+
 @pytest.mark.parametrize("shots", [1000])
 def test_probability(token, tol, shots):
     """Test that the probs function works."""
diff --git a/tests/test_integration.py b/tests/test_integration.py
index 88dbacb18..e45ee2762 100644
--- a/tests/test_integration.py
+++ b/tests/test_integration.py
@@ -8,6 +8,7 @@
 import qiskit.providers.aer as aer
 
 from pennylane_qiskit import AerDevice, BasicAerDevice
+from pennylane_qiskit.qiskit_device import QiskitDevice
 
 from conftest import state_backends
 
@@ -54,6 +55,10 @@ def test_args(self):
     @pytest.mark.parametrize("shots", [None, 8192])
     def test_one_qubit_circuit(self, shots, d, backend, tol):
         """Test that devices provide correct result for a simple circuit"""
+        if (d[0] == "qiskit.aer" and "aer" not in backend) \
+          or (d[0] == "qiskit.basicaer" and "aer" in backend):
+            pytest.skip("Only the AerSimulator is supported on AerDevice")
+
         if backend not in state_backends and shots is None:
             pytest.skip("Hardware simulators do not support analytic mode")
 
@@ -75,7 +80,7 @@ def circuit(x, y, z):
 
     @pytest.mark.parametrize("d", pldevices)
     @pytest.mark.parametrize("shots", [8192])
-    def test_one_qubit_circuit(self, shots, d, backend, tol):
+    def test_basis_state_and_rot(self, shots, d, backend, tol):
         """Integration test for the BasisState and Rot operations for non-analytic mode."""
 
         if (d[0] == "qiskit.aer" and "aer" not in backend) \
@@ -479,3 +484,82 @@ def circuit():
             return qml.expval(qml.PauliZ(wires=0))
 
         assert circuit() == -1
+
+class TestBatchExecution:
+    """Test the devices work correctly with the batch execution pipeline."""
+
+    @pytest.mark.parametrize("d", pldevices)
+    @pytest.mark.parametrize("shots", [None, 8192])
+    def test_one_qubit_circuit_batch_params(self, shots, d, backend, tol, mocker):
+        """Test that devices provide correct result for a simple circuit using
+        the batch_params transform."""
+        if (d[0] == "qiskit.aer" and "aer" not in backend) \
+          or (d[0] == "qiskit.basicaer" and "aer" in backend):
+            pytest.skip("Only the AerSimulator is supported on AerDevice")
+
+        if backend not in state_backends and shots is None:
+            pytest.skip("Hardware simulators do not support analytic mode")
+
+        dev = qml.device(d[0], wires=1, backend=backend, shots=shots)
+
+        # Batch the input parameters
+        batch_dim = 3
+        a = np.linspace(0, 0.543, batch_dim)
+        b = np.linspace(0, 0.123, batch_dim)
+        c = np.linspace(0, 0.987, batch_dim)
+
+        spy1 = mocker.spy(QiskitDevice, "batch_execute")
+        spy2 = mocker.spy(dev.backend, "run")
+
+        @qml.batch_params
+        @qml.qnode(dev)
+        def circuit(x, y, z):
+            """Reference QNode"""
+            qml.PauliX(0)
+            qml.Hadamard(wires=0)
+            qml.Rot(x, y, z, wires=0)
+            return qml.expval(qml.PauliZ(0))
+
+        assert np.allclose(circuit(a, b, c), np.cos(a) * np.sin(b), **tol)
+
+        # Check that QiskitDevice.batch_execute was called
+        assert spy1.call_count == 1
+        assert spy2.call_count == 1
+
+    @pytest.mark.parametrize("d", pldevices)
+    @pytest.mark.parametrize("shots", [None, 8192])
+    def test_batch_execute_parameter_shift(self, shots, d, backend, tol, mocker):
+        """Test that devices provide correct result computing the gradient of a
+        circuit using the parameter-shift rule and the batch execution pipeline."""
+        if (d[0] == "qiskit.aer" and "aer" not in backend) \
+          or (d[0] == "qiskit.basicaer" and "aer" in backend):
+            pytest.skip("Only the AerSimulator is supported on AerDevice")
+
+        if backend not in state_backends and shots is None:
+            pytest.skip("Hardware simulators do not support analytic mode")
+
+        dev = qml.device(d[0], wires=3, backend=backend, shots=shots)
+
+        spy1 = mocker.spy(QiskitDevice, "batch_execute")
+        spy2 = mocker.spy(dev.backend, "run")
+
+        @qml.qnode(dev, diff_method="parameter-shift")
+        def circuit(x, y):
+            qml.RX(x, wires=[0])
+            qml.RY(y, wires=[1])
+            qml.CNOT(wires=[0, 1])
+            return qml.expval(qml.PauliZ(0) @ qml.PauliX(1) @ qml.PauliZ(2))
+
+        x = qml.numpy.array(0.543, requires_grad=True)
+        y = qml.numpy.array(0.123, requires_grad=True)
+
+        res = qml.grad(circuit)(x,y)
+        expected = np.array([[-np.sin(y) * np.sin(x), np.cos(y) * np.cos(x)]])
+        assert np.allclose(res, expected, **tol)
+
+        # Check that QiskitDevice.batch_execute was called once
+        assert spy1.call_count == 1
+
+        # Check that run was called twice: for the partial derivatives and for
+        # running the circuit
+        assert spy2.call_count == 2
diff --git a/tests/test_qiskit_device.py b/tests/test_qiskit_device.py
index 59862182a..ad987ca88 100644
--- a/tests/test_qiskit_device.py
+++ b/tests/test_qiskit_device.py
@@ -3,6 +3,7 @@
 
 import pennylane as qml
 from pennylane_qiskit import AerDevice
+from pennylane_qiskit.qiskit_device import QiskitDevice
 import qiskit.providers.aer.noise as noise
 
 test_transpile_options = [
@@ -98,3 +99,81 @@ def test_backend_options_cleaned(self):
 
         dev2 = qml.device("qiskit.aer", wires=2)
         assert dev2.backend.options.get("noise_model") is None
+
+@pytest.mark.parametrize("shots", [None])
+class TestBatchExecution:
+    """Tests for the batch_execute method."""
+
+    with qml.tape.QuantumTape() as tape1:
+        qml.PauliX(wires=0)
+        qml.expval(qml.PauliZ(wires=0)), qml.expval(qml.PauliZ(wires=1))
+
+    with qml.tape.JacobianTape() as tape2:
+        qml.PauliX(wires=0)
+        qml.expval(qml.PauliZ(wires=0))
+
+    @pytest.mark.parametrize("n_tapes", [1, 2, 3])
+    def test_calls_to_execute(self, device, n_tapes, mocker):
+        """Tests that only the device's dedicated batch execute method is
+        called and not the general execute method."""
+
+        dev = device(2)
+        spy = mocker.spy(QiskitDevice, "execute")
+
+        tapes = [self.tape1] * n_tapes
+        dev.batch_execute(tapes)
+
+        # Check that QiskitDevice.execute was not called
+        assert spy.call_count == 0
+
+    @pytest.mark.parametrize("n_tapes", [1, 2, 3])
+    def test_calls_to_reset(self, n_tapes, mocker, device):
+        """Tests that the device's reset method is called the correct number of
+        times."""
+
+        dev = device(2)
+        spy = mocker.spy(QiskitDevice, "reset")
+
+        tapes = [self.tape1] * n_tapes
+        dev.batch_execute(tapes)
+
+        assert spy.call_count == n_tapes
+
+    def test_result(self, device, tol):
+        """Tests that the result has the correct shape and entry types."""
+        dev = device(2)
+        tapes = [self.tape1, self.tape2]
+        res = dev.batch_execute(tapes)
+
+        # We're calling device methods directly, need to reset before the next
+        # execution
+        dev.reset()
+        tape1_expected = dev.execute(self.tape1)
+
+        dev.reset()
+        tape2_expected = dev.execute(self.tape2)
+
+        assert len(res) == 2
+        assert np.allclose(
+            res[0], tape1_expected, atol=0
+        )
+
+        assert np.allclose(
+            res[1], tape2_expected, atol=0
+        )
+
+    def test_result_empty_tape(self, device, tol):
+        """Tests that the result has the correct shape and entry types for empty tapes."""
+        dev = device(2)
+
+        empty_tape = qml.tape.QuantumTape()
+        tapes = [empty_tape] * 3
+        res = dev.batch_execute(tapes)
+
+        # We're calling device methods directly, need to reset before the next
+        # execution
+        dev.reset()
+        assert len(res) == 3
+        assert np.allclose(
+            res[0], dev.execute(empty_tape), atol=0
+        )