diff --git a/qiskit_experiments/library/randomized_benchmarking/clifford_synthesis.py b/qiskit_experiments/library/randomized_benchmarking/clifford_synthesis.py
new file mode 100644
index 0000000000..29873fb37f
--- /dev/null
+++ b/qiskit_experiments/library/randomized_benchmarking/clifford_synthesis.py
@@ -0,0 +1,101 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Clifford synthesis plugins for randomized benchmarking
+"""
+from __future__ import annotations
+
+from typing import Sequence
+
+from qiskit.circuit import QuantumCircuit, Operation
+from qiskit.circuit.equivalence_library import SessionEquivalenceLibrary as sel
+from qiskit.exceptions import QiskitError
+from qiskit.synthesis.clifford import synth_clifford_full
+from qiskit.transpiler import PassManager, CouplingMap, Layout, Target
+from qiskit.transpiler.passes import (
+    SabreSwap,
+    LayoutTransformation,
+    BasisTranslator,
+    CheckGateDirection,
+    GateDirection,
+    Optimize1qGatesDecomposition,
+)
+from qiskit.transpiler.passes.synthesis.plugin import HighLevelSynthesisPlugin
+
+
+class RBDefaultCliffordSynthesis(HighLevelSynthesisPlugin):
+    """Default Clifford synthesis plugin for randomized benchmarking."""
+
+    def run(
+        self,
+        high_level_object: Operation,
+        coupling_map: CouplingMap | None = None,
+        target: Target | None = None,
+        qubits: Sequence | None = None,
+        **options,
+    ) -> QuantumCircuit:
+        """Run synthesis for the given Clifford.
+
+        Args:
+            high_level_object: The operation to synthesize to a
+                :class:`~qiskit.circuit.QuantumCircuit` object.
+            coupling_map: The reduced coupling map of the backend. For example,
+                if physical qubits [5, 6, 7] to be benchmarked is connected
+                as 5 - 7 - 6 linearly, the reduced coupling map is 0 - 2 - 1.
+            target: A target representing the target backend, which will be ignored in this plugin.
+            qubits: List of physical qubits over which the operation is defined,
+                which will be ignored in this plugin.
+            options: Additional method-specific optional kwargs,
+                which must include ``basis_gates``, basis gates to be used for the synthesis.
+
+        Returns:
+            The quantum circuit representation of the Operation
+                when successful, and ``None`` otherwise.
+
+        Raises:
+            QiskitError: If basis_gates is not supplied.
+        """
+        # synthesize cliffords
+        circ = synth_clifford_full(high_level_object)
+
+        # post processing to comply with basis gates and coupling map
+        if coupling_map is None:  # Sabre does not work with coupling_map=None
+            return circ
+
+        basis_gates = options.get("basis_gates", None)
+        if basis_gates is None:
+            raise QiskitError("basis_gates are required to run this synthesis plugin")
+
+        basis_gates = list(basis_gates)
+
+        # Run Sabre routing and undo the layout change
+        # assuming Sabre routing does not change the initial layout.
+        # And then decompose swap gates, fix 2q-gate direction and optimize 1q gates
+        initial_layout = Layout.generate_trivial_layout(*circ.qubits)
+        undo_layout_change = LayoutTransformation(
+            coupling_map=coupling_map, from_layout="final_layout", to_layout=initial_layout
+        )
+
+        def _direction_condition(property_set):
+            return not property_set["is_direction_mapped"]
+
+        pm = PassManager(
+            [
+                SabreSwap(coupling_map),
+                undo_layout_change,
+                BasisTranslator(sel, basis_gates),
+                CheckGateDirection(coupling_map),
+            ]
+        )
+        pm.append([GateDirection(coupling_map)], condition=_direction_condition)
+        pm.append([Optimize1qGatesDecomposition(basis=basis_gates)])
+        return pm.run(circ)
diff --git a/qiskit_experiments/library/randomized_benchmarking/clifford_utils.py b/qiskit_experiments/library/randomized_benchmarking/clifford_utils.py
index f6ab757b5a..a081db7006 100644
--- a/qiskit_experiments/library/randomized_benchmarking/clifford_utils.py
+++ b/qiskit_experiments/library/randomized_benchmarking/clifford_utils.py
@@ -29,8 +29,12 @@
 from qiskit.compiler import transpile
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info import Clifford, random_clifford
+from qiskit.transpiler import CouplingMap, PassManager
+from qiskit.transpiler.passes.synthesis.high_level_synthesis import HLSConfig, HighLevelSynthesis
 from qiskit.utils.deprecation import deprecate_func
 
+DEFAULT_SYNTHESIS_METHOD = "rb_default"
+
 _DATA_FOLDER = os.path.join(os.path.dirname(__file__), "data")
 
 _CLIFFORD_COMPOSE_1Q = np.load(f"{_DATA_FOLDER}/clifford_compose_1q.npz")["table"]
@@ -110,37 +114,141 @@ def _circuit_compose(
     return self
 
 
-def _truncate_inactive_qubits(
-    circ: QuantumCircuit, active_qubits: Sequence[Qubit]
+def _synthesize_clifford(
+    clifford: Clifford,
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]] = None,
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> QuantumCircuit:
-    res = QuantumCircuit(active_qubits, name=circ.name, metadata=circ.metadata)
-    for inst in circ:
-        if all(q in active_qubits for q in inst.qubits):
-            res.append(inst)
-    res.calibrations = circ.calibrations
-    return res
+    """Synthesize a circuit of a Clifford element. The resulting circuit contains only
+    ``basis_gates`` and it complies with ``coupling_tuple``.
+
+    Args:
+        clifford: Clifford element to be converted
+        basis_gates: basis gates to use in the conversion
+        coupling_tuple: coupling map to use in the conversion in the form of tuple of edges
+        synthesis_method: conversion algorithm name
+
+    Returns:
+        Synthesized circuit
+    """
+    qc = QuantumCircuit(clifford.num_qubits, name=str(clifford))
+    qc.append(clifford, qc.qubits)
+    return _synthesize_clifford_circuit(
+        qc,
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    )
 
 
 def _synthesize_clifford_circuit(
-    circuit: QuantumCircuit, basis_gates: Tuple[str]
+    circuit: QuantumCircuit,
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]] = None,
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> QuantumCircuit:
-    # synthesizes clifford circuits using given basis gates, for use during
-    # custom transpilation during RB circuit generation.
-    return transpile(circuit, basis_gates=list(basis_gates), optimization_level=1)
+    """Convert a Clifford circuit into one composed of ``basis_gates`` with
+    satisfying ``coupling_tuple`` using the specified synthesis method.
+
+    Args:
+        circuit: Clifford circuit to be converted
+        basis_gates: basis gates to use in the conversion
+        coupling_tuple: coupling map to use in the conversion in the form of tuple of edges
+        synthesis_method: name of Clifford synthesis algorithm to use
+
+    Returns:
+        Synthesized circuit
+    """
+    if basis_gates:
+        basis_gates = list(basis_gates)
+    coupling_map = CouplingMap(coupling_tuple) if coupling_tuple else None
+
+    # special handling for 1q or 2q case for speed
+    if circuit.num_qubits <= 2:
+        if synthesis_method == DEFAULT_SYNTHESIS_METHOD:
+            return transpile(
+                circuit,
+                basis_gates=basis_gates,
+                coupling_map=coupling_map,
+                optimization_level=1,
+            )
+        else:
+            # Provided custom synthesis method, re-synthesize Clifford circuit
+            # convert the circuit back to a Clifford object and then call the synthesis plugin
+            new_circuit = QuantumCircuit(circuit.num_qubits, name=circuit.name)
+            new_circuit.append(Clifford(circuit), new_circuit.qubits)
+            circuit = new_circuit
+
+    # for 3q+ or custom synthesis method, synthesizes clifford circuit
+    hls_config = HLSConfig(clifford=[(synthesis_method, {"basis_gates": basis_gates})])
+    pm = PassManager([HighLevelSynthesis(hls_config=hls_config, coupling_map=coupling_map)])
+    circuit = pm.run(circuit)
+    return circuit
 
 
-@lru_cache(maxsize=None)
+@lru_cache(maxsize=256)
 def _clifford_1q_int_to_instruction(
-    num: Integral, basis_gates: Optional[Tuple[str]]
+    num: Integral,
+    basis_gates: Optional[Tuple[str]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> Instruction:
-    return CliffordUtils.clifford_1_qubit_circuit(num, basis_gates).to_instruction()
+    return CliffordUtils.clifford_1_qubit_circuit(
+        num, basis_gates=basis_gates, synthesis_method=synthesis_method
+    ).to_instruction()
 
 
 @lru_cache(maxsize=11520)
 def _clifford_2q_int_to_instruction(
-    num: Integral, basis_gates: Optional[Tuple[str]]
+    num: Integral,
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> Instruction:
-    return CliffordUtils.clifford_2_qubit_circuit(num, basis_gates).to_instruction()
+    return CliffordUtils.clifford_2_qubit_circuit(
+        num,
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    ).to_instruction()
+
+
+def _hash_cliff(cliff):
+    return cliff.tableau.tobytes(), cliff.tableau.shape
+
+
+def _dehash_cliff(cliff_hash):
+    tableau = np.frombuffer(cliff_hash[0], dtype=bool).reshape(cliff_hash[1])
+    return Clifford(tableau)
+
+
+def _clifford_to_instruction(
+    clifford: Clifford,
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
+) -> Instruction:
+    return _cached_clifford_to_instruction(
+        _hash_cliff(clifford),
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    )
+
+
+@lru_cache(maxsize=256)
+def _cached_clifford_to_instruction(
+    cliff_hash: Tuple[str, Tuple[int, int]],
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
+) -> Instruction:
+    return _synthesize_clifford(
+        _dehash_cliff(cliff_hash),
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    ).to_instruction()
 
 
 # The classes VGate and WGate are not actually used in the code - we leave them here to give
@@ -254,7 +362,12 @@ def random_clifford_circuits(
 
     @classmethod
     @lru_cache(maxsize=24)
-    def clifford_1_qubit_circuit(cls, num, basis_gates: Optional[Tuple[str, ...]] = None):
+    def clifford_1_qubit_circuit(
+        cls,
+        num,
+        basis_gates: Optional[Tuple[str, ...]] = None,
+        synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
+    ):
         """Return the 1-qubit clifford circuit corresponding to ``num``,
         where ``num`` is between 0 and 23.
         """
@@ -275,20 +388,28 @@ def clifford_1_qubit_circuit(cls, num, basis_gates: Optional[Tuple[str, ...]] =
             qc.z(0)
 
         if basis_gates:
-            qc = _synthesize_clifford_circuit(qc, basis_gates)
+            qc = _synthesize_clifford_circuit(qc, basis_gates, synthesis_method=synthesis_method)
 
         return qc
 
     @classmethod
     @lru_cache(maxsize=11520)
-    def clifford_2_qubit_circuit(cls, num, basis_gates: Optional[Tuple[str, ...]] = None):
+    def clifford_2_qubit_circuit(
+        cls,
+        num,
+        basis_gates: Optional[Tuple[str, ...]] = None,
+        coupling_tuple: Optional[Tuple[Tuple[int, int]]] = None,
+        synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
+    ):
         """Return the 2-qubit clifford circuit corresponding to `num`
         where `num` is between 0 and 11519.
         """
         qc = QuantumCircuit(2, name=f"Clifford-2Q({num})")
         for layer, idx in enumerate(_layer_indices_from_num(num)):
             if basis_gates:
-                layer_circ = _transformed_clifford_layer(layer, idx, basis_gates)
+                layer_circ = _transformed_clifford_layer(
+                    layer, idx, basis_gates, coupling_tuple, synthesis_method=synthesis_method
+                )
             else:
                 layer_circ = _CLIFFORD_LAYER[layer][idx]
             _circuit_compose(qc, layer_circ, qubits=(0, 1))
@@ -578,13 +699,22 @@ def _clifford_2q_nums_from_2q_circuit(qc: QuantumCircuit) -> Iterable[Integral]:
 ]
 
 
-@lru_cache(maxsize=None)
+@lru_cache(maxsize=256)
 def _transformed_clifford_layer(
-    layer: int, index: Integral, basis_gates: Tuple[str, ...]
+    layer: int,
+    index: Integral,
+    basis_gates: Tuple[str, ...],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> QuantumCircuit:
     # Return the index-th quantum circuit of the layer translated with the basis_gates.
     # The result is cached for speed.
-    return _synthesize_clifford_circuit(_CLIFFORD_LAYER[layer][index], basis_gates)
+    return _synthesize_clifford_circuit(
+        _CLIFFORD_LAYER[layer][index],
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    )
 
 
 def _num_from_layer_indices(triplet: Tuple[Integral, Integral, Integral]) -> Integral:
diff --git a/qiskit_experiments/library/randomized_benchmarking/interleaved_rb_experiment.py b/qiskit_experiments/library/randomized_benchmarking/interleaved_rb_experiment.py
index 0d105ece0f..a1b2c8d179 100644
--- a/qiskit_experiments/library/randomized_benchmarking/interleaved_rb_experiment.py
+++ b/qiskit_experiments/library/randomized_benchmarking/interleaved_rb_experiment.py
@@ -14,20 +14,18 @@
 """
 import itertools
 import warnings
-from typing import Union, Iterable, Optional, List, Sequence, Tuple
+from typing import Union, Iterable, Optional, List, Sequence, Dict, Any
 
 from numpy.random import Generator
 from numpy.random.bit_generator import BitGenerator, SeedSequence
 
 from qiskit.circuit import QuantumCircuit, Instruction, Gate, Delay
-from qiskit.compiler import transpile
 from qiskit.exceptions import QiskitError
 from qiskit.providers.backend import Backend
 from qiskit.quantum_info import Clifford
-from qiskit.transpiler.exceptions import TranspilerError
 from qiskit_experiments.framework import Options
 from qiskit_experiments.framework.backend_timing import BackendTiming
-from .clifford_utils import _truncate_inactive_qubits
+from .clifford_utils import _synthesize_clifford
 from .clifford_utils import num_from_1q_circuit, num_from_2q_circuit
 from .interleaved_rb_analysis import InterleavedRBAnalysis
 from .standard_rb import StandardRB, SequenceElementType
@@ -73,8 +71,7 @@ def __init__(
         Args:
             interleaved_element: The element to interleave,
                     given either as a Clifford element, gate, delay or circuit.
-                    If the element contains any non-basis gates,
-                    it will be transpiled with ``transpiled_options`` of this experiment.
+                    All instructions in the element must be supported in the ``backend``(``target``).
                     If it is/contains a delay, its duration and unit must comply with
                     the timing constraints of the ``backend``
                     (:class:`~qiskit_experiments.framework.backend_timing.BackendTiming`
@@ -186,8 +183,7 @@ def circuits(self) -> List[QuantumCircuit]:
             A list of :class:`QuantumCircuit`.
 
         Raises:
-            QiskitError: If the ``interleaved_element`` provided to the constructor
-                cannot be transpiled.
+            QiskitError: If interleaved_element has non-supported instruction in the backend.
         """
         # Convert interleaved element to transpiled circuit operation and store it for speed
         self.__set_up_interleaved_op()
@@ -223,20 +219,22 @@ def circuits(self) -> List[QuantumCircuit]:
         return list(itertools.chain.from_iterable(zip(reference_circuits, interleaved_circuits)))
 
     def _to_instruction(
-        self, elem: SequenceElementType, basis_gates: Optional[Tuple[str]] = None
+        self,
+        elem: SequenceElementType,
+        synthesis_options: Dict[str, Optional[Any]],
     ) -> Instruction:
         if elem is self._interleaved_cliff:
             return self._interleaved_op
 
-        return super()._to_instruction(elem, basis_gates)
+        return super()._to_instruction(elem, synthesis_options)
 
     def __set_up_interleaved_op(self) -> None:
         # Convert interleaved element to transpiled circuit operation and store it for speed
         self._interleaved_op = self._interleaved_element
-        basis_gates = self._get_basis_gates()
         # Convert interleaved element to circuit
         if isinstance(self._interleaved_op, Clifford):
-            self._interleaved_op = self._interleaved_op.to_circuit()
+            opts = self._get_synthesis_options()
+            self._interleaved_op = _synthesize_clifford(self._interleaved_op, **opts)
 
         if isinstance(self._interleaved_op, QuantumCircuit):
             interleaved_circ = self._interleaved_op
@@ -246,22 +244,17 @@ def __set_up_interleaved_op(self) -> None:
         else:  # Delay
             interleaved_circ = []
 
-        if basis_gates and any(i.operation.name not in basis_gates for i in interleaved_circ):
-            # Transpile circuit with non-basis gates and remove idling qubits
-            try:
-                interleaved_circ = transpile(
-                    interleaved_circ, self.backend, **vars(self.transpile_options)
+        # Validate if all instructions in the interleaved circuit are supported in the backend
+        if self.backend and hasattr(self.backend, "target"):
+            for inst in interleaved_circ:
+                qargs = tuple(
+                    self.physical_qubits[interleaved_circ.find_bit(q).index] for q in inst.qubits
                 )
-            except TranspilerError as err:
-                raise QiskitError("Failed to transpile interleaved_element.") from err
-            interleaved_circ = _truncate_inactive_qubits(
-                interleaved_circ, active_qubits=interleaved_circ.qubits[: self.num_qubits]
-            )
-            # Convert transpiled circuit to operation
-            if len(interleaved_circ) == 1:
-                self._interleaved_op = interleaved_circ.data[0].operation
-            else:
-                self._interleaved_op = interleaved_circ
+                if not self.backend.target.instruction_supported(inst.operation.name, qargs):
+                    raise QiskitError(
+                        f"{inst.operation.name} in interleaved element is not supported"
+                        f" on qubits {qargs} in the backend."
+                    )
 
         # Store interleaved operation as Instruction
         if isinstance(self._interleaved_op, QuantumCircuit):
diff --git a/qiskit_experiments/library/randomized_benchmarking/standard_rb.py b/qiskit_experiments/library/randomized_benchmarking/standard_rb.py
index da03d37bf8..9a7795155a 100644
--- a/qiskit_experiments/library/randomized_benchmarking/standard_rb.py
+++ b/qiskit_experiments/library/randomized_benchmarking/standard_rb.py
@@ -12,17 +12,18 @@
 """
 Standard RB Experiment class.
 """
-import logging
 import functools
+import logging
 from collections import defaultdict
 from numbers import Integral
-from typing import Union, Iterable, Optional, List, Sequence, Tuple
+from typing import Union, Iterable, Optional, List, Sequence, Dict, Any
 
 import numpy as np
+import rustworkx as rx
 from numpy.random import Generator, default_rng
 from numpy.random.bit_generator import BitGenerator, SeedSequence
 
-from qiskit.circuit import CircuitInstruction, QuantumCircuit, Instruction, Barrier
+from qiskit.circuit import CircuitInstruction, QuantumCircuit, Instruction, Barrier, Gate
 from qiskit.exceptions import QiskitError
 from qiskit.providers import BackendV2Converter
 from qiskit.providers.backend import Backend, BackendV1, BackendV2
@@ -30,18 +31,18 @@
 from qiskit.quantum_info import Clifford
 from qiskit.quantum_info.random import random_clifford
 from qiskit.transpiler import CouplingMap
-
 from qiskit_experiments.framework import BaseExperiment, Options
 from qiskit_experiments.framework.restless_mixin import RestlessMixin
-
 from .clifford_utils import (
     CliffordUtils,
+    DEFAULT_SYNTHESIS_METHOD,
     compose_1q,
     compose_2q,
     inverse_1q,
     inverse_2q,
     _clifford_1q_int_to_instruction,
     _clifford_2q_int_to_instruction,
+    _clifford_to_instruction,
     _transpile_clifford_circuit,
 )
 from .rb_analysis import RBAnalysis
@@ -147,6 +148,8 @@ def _default_experiment_options(cls) -> Options:
             full_sampling (bool): If True all Cliffords are independently sampled for
                 all lengths. If False for sample of lengths longer sequences are constructed
                 by appending additional Clifford samples to shorter sequences.
+            clifford_synthesis_method (str): The name of the Clifford synthesis plugin to use
+                for building circuits of RB sequences.
         """
         options = super()._default_experiment_options()
         options.update_options(
@@ -154,6 +157,7 @@ def _default_experiment_options(cls) -> Options:
             num_samples=None,
             seed=None,
             full_sampling=None,
+            clifford_synthesis_method=DEFAULT_SYNTHESIS_METHOD,
         )
 
         return options
@@ -210,64 +214,58 @@ def _sample_sequences(self) -> List[Sequence[SequenceElementType]]:
 
         return sequences
 
-    def _get_basis_gates(self) -> Optional[Tuple[str, ...]]:
-        """Get sorted basis gates to use in basis transformation during circuit generation.
-
-        - Return None if this experiment is an RB with 3 or more qubits.
-        - Return None if no basis gates are supplied via ``backend`` or ``transpile_options``.
-        - Return None if all 2q-gates supported on the physical qubits of the backend are one-way
-        directed (e.g. cx(0, 1) is supported but cx(1, 0) is not supported).
+    def _get_synthesis_options(self) -> Dict[str, Optional[Any]]:
+        """Get options for Clifford synthesis from the backend information as a dictionary.
 
-        In all those case when None are returned, basis transformation will be skipped in the
-        circuit generation step (i.e. :meth:`circuits`) and it will be done in the successive
-        transpilation step (i.e. :meth:`_transpiled_circuits`) that calls :func:`transpile`.
+        The options include:
+        - "basis_gates": Sorted basis gate names.
+            Return None if no basis gates are supplied via ``backend`` or ``transpile_options``.
+        - "coupling_tuple": Reduced coupling map in the form of tuple of edges in the coupling graph.
+            Return None if no coupling map are supplied via ``backend`` or ``transpile_options``.
 
         Returns:
-            Sorted basis gate names.
+            Synthesis options as a dictionary.
         """
-        # 3 or more qubits case: Return None (skip basis transformation in circuit generation)
-        if self.num_qubits > 2:
-            return None
-
-        # 1 qubit case: Return all basis gates (or None if no basis gates are supplied)
-        if self.num_qubits == 1:
-            basis_gates = self.transpile_options.get("basis_gates", None)
-            if not basis_gates and self.backend:
-                if isinstance(self.backend, BackendV2):
-                    basis_gates = self.backend.operation_names
-                elif isinstance(self.backend, BackendV1):
-                    basis_gates = self.backend.configuration().basis_gates
-            return tuple(sorted(basis_gates)) if basis_gates else None
-
-        def is_bidirectional(coupling_map):
-            if coupling_map is None:
-                # None for a coupling map implies all-to-all coupling
-                return True
-            return len(coupling_map.reduce(self.physical_qubits).get_edges()) == 2
-
-        # 2 qubits case: Return all basis gates except for one-way directed 2q-gates.
-        # Return None if there is no bidirectional 2q-gates in basis gates.
-        if self.num_qubits == 2:
-            basis_gates = self.transpile_options.get("basis_gates", [])
-            if not basis_gates and self.backend:
-                if isinstance(self.backend, BackendV2) and self.backend.target:
-                    has_bidirectional_2q_gates = False
-                    for op_name in self.backend.target:
-                        if self.backend.target.operation_from_name(op_name).num_qubits == 2:
-                            if is_bidirectional(self.backend.target.build_coupling_map(op_name)):
-                                has_bidirectional_2q_gates = True
-                            else:
-                                continue
-                        basis_gates.append(op_name)
-                    if not has_bidirectional_2q_gates:
-                        basis_gates = None
-                elif isinstance(self.backend, BackendV1):
-                    cmap = self.backend.configuration().coupling_map
-                    if cmap is None or is_bidirectional(CouplingMap(cmap)):
-                        basis_gates = self.backend.configuration().basis_gates
-            return tuple(sorted(basis_gates)) if basis_gates else None
-
-        return None
+        basis_gates = self.transpile_options.get("basis_gates", [])
+        coupling_map = self.transpile_options.get("coupling_map", None)
+        if coupling_map:
+            coupling_map = coupling_map.reduce(self.physical_qubits)
+        if not (basis_gates and coupling_map) and self.backend:
+            if isinstance(self.backend, BackendV2) and "simulator" in self.backend.name:
+                basis_gates = basis_gates if basis_gates else self.backend.target.operation_names
+                coupling_map = coupling_map if coupling_map else None
+            elif isinstance(self.backend, BackendV2):
+                gate_ops = [op for op in self.backend.target.operations if isinstance(op, Gate)]
+                backend_basis_gates = [op.name for op in gate_ops if op.num_qubits != 2]
+                backend_cmap = None
+                for op in gate_ops:
+                    if op.num_qubits != 2:
+                        continue
+                    cmap = self.backend.target.build_coupling_map(op.name)
+                    if cmap is None:
+                        backend_basis_gates.append(op.name)
+                    else:
+                        reduced = cmap.reduce(self.physical_qubits)
+                        if rx.is_weakly_connected(reduced.graph):
+                            backend_basis_gates.append(op.name)
+                            backend_cmap = reduced
+                            # take the first non-global 2q gate if backend has multiple 2q gates
+                            break
+                basis_gates = basis_gates if basis_gates else backend_basis_gates
+                coupling_map = coupling_map if coupling_map else backend_cmap
+            elif isinstance(self.backend, BackendV1):
+                backend_basis_gates = self.backend.configuration().basis_gates
+                backend_cmap = self.backend.configuration().coupling_map
+                if backend_cmap:
+                    backend_cmap = CouplingMap(backend_cmap).reduce(self.physical_qubits)
+                basis_gates = basis_gates if basis_gates else backend_basis_gates
+                coupling_map = coupling_map if coupling_map else backend_cmap
+
+        return {
+            "basis_gates": tuple(sorted(basis_gates)) if basis_gates else None,
+            "coupling_tuple": tuple(sorted(coupling_map.get_edges())) if coupling_map else None,
+            "synthesis_method": self.experiment_options["clifford_synthesis_method"],
+        }
 
     def _sequences_to_circuits(
         self, sequences: List[Sequence[SequenceElementType]]
@@ -277,7 +275,7 @@ def _sequences_to_circuits(
         Returns:
             A list of RB circuits.
         """
-        basis_gates = self._get_basis_gates()
+        synthesis_opts = self._get_synthesis_options()
         # Circuit generation
         circuits = []
         for i, seq in enumerate(sequences):
@@ -289,7 +287,7 @@ def _sequences_to_circuits(
 
             circ = QuantumCircuit(self.num_qubits)
             for elem in seq:
-                circ.append(self._to_instruction(elem, basis_gates), circ.qubits)
+                circ.append(self._to_instruction(elem, synthesis_opts), circ.qubits)
                 circ._append(CircuitInstruction(Barrier(self.num_qubits), circ.qubits))
 
             # Compute inverse, compute only the difference from the previous shorter sequence
@@ -297,7 +295,7 @@ def _sequences_to_circuits(
             prev_seq = seq
             inv = self.__adjoint_clifford(prev_elem)
 
-            circ.append(self._to_instruction(inv, basis_gates), circ.qubits)
+            circ.append(self._to_instruction(inv, synthesis_opts), circ.qubits)
             circ.measure_all()  # includes insertion of the barrier before measurement
             circuits.append(circ)
         return circuits
@@ -309,20 +307,38 @@ def __sample_sequence(self, length: int, rng: Generator) -> Sequence[SequenceEle
             return rng.integers(CliffordUtils.NUM_CLIFFORD_1_QUBIT, size=length)
         if self.num_qubits == 2:
             return rng.integers(CliffordUtils.NUM_CLIFFORD_2_QUBIT, size=length)
-        # Return circuit object instead of Clifford object for 3 or more qubits case for speed
-        return [random_clifford(self.num_qubits, rng).to_circuit() for _ in range(length)]
+        # Return Clifford object for 3 or more qubits case
+        return [random_clifford(self.num_qubits, rng) for _ in range(length)]
 
     def _to_instruction(
-        self, elem: SequenceElementType, basis_gates: Optional[Tuple[str, ...]] = None
+        self,
+        elem: SequenceElementType,
+        synthesis_options: Dict[str, Optional[Any]],
     ) -> Instruction:
+        """Return the instruction of a Clifford element.
+
+        The resulting instruction contains a circuit definition with ``basis_gates`` and
+        it complies with ``coupling_tuple``, which is specified in ``synthesis_options``.
+
+        Args:
+            elem: a Clifford element to be converted
+            synthesis_options: options for synthesizing the Clifford element
+
+        Returns:
+            Converted instruction
+        """
         # Switching for speed up
         if isinstance(elem, Integral):
             if self.num_qubits == 1:
-                return _clifford_1q_int_to_instruction(elem, basis_gates)
+                return _clifford_1q_int_to_instruction(
+                    elem,
+                    basis_gates=synthesis_options["basis_gates"],
+                    synthesis_method=synthesis_options["synthesis_method"],
+                )
             if self.num_qubits == 2:
-                return _clifford_2q_int_to_instruction(elem, basis_gates)
+                return _clifford_2q_int_to_instruction(elem, **synthesis_options)
 
-        return elem.to_instruction()
+        return _clifford_to_instruction(elem, **synthesis_options)
 
     def __identity_clifford(self) -> SequenceElementType:
         if self.num_qubits <= 2:
@@ -336,11 +352,11 @@ def __compose_clifford_seq(
             return functools.reduce(
                 compose_1q if self.num_qubits == 1 else compose_2q, elements, base_elem
             )
-        # 3 or more qubits: compose Clifford from circuits for speed
-        circ = QuantumCircuit(self.num_qubits)
+        # 3 or more qubits
+        res = base_elem
         for elem in elements:
-            circ.compose(elem, inplace=True)
-        return base_elem.compose(Clifford.from_circuit(circ))
+            res = res.compose(elem)
+        return res
 
     def __adjoint_clifford(self, op: SequenceElementType) -> SequenceElementType:
         if self.num_qubits == 1:
@@ -354,20 +370,21 @@ def __adjoint_clifford(self, op: SequenceElementType) -> SequenceElementType:
     def _transpiled_circuits(self) -> List[QuantumCircuit]:
         """Return a list of experiment circuits, transpiled."""
         has_custom_transpile_option = (
-            not set(vars(self.transpile_options)).issubset({"basis_gates", "optimization_level"})
+            not set(vars(self.transpile_options)).issubset(
+                {"basis_gates", "coupling_map", "optimization_level"}
+            )
             or self.transpile_options.get("optimization_level", 1) != 1
         )
-        has_no_undirected_2q_basis = self._get_basis_gates() is None
-        if self.num_qubits > 2 or has_custom_transpile_option or has_no_undirected_2q_basis:
+        if has_custom_transpile_option:
             transpiled = super()._transpiled_circuits()
         else:
             transpiled = [
                 _transpile_clifford_circuit(circ, physical_qubits=self.physical_qubits)
                 for circ in self.circuits()
             ]
-            # Set custom calibrations provided in backend
-            if isinstance(self.backend, BackendV2):
-                qargs_patterns = [self.physical_qubits]  # for self.num_qubits == 1
+            # Set custom calibrations provided in backend (excluding simulators)
+            if isinstance(self.backend, BackendV2) and "simulator" not in self.backend.name:
+                qargs_patterns = [self.physical_qubits]  # for 1q or 3q+ case
                 if self.num_qubits == 2:
                     qargs_patterns = [
                         (self.physical_qubits[0],),
@@ -376,10 +393,12 @@ def _transpiled_circuits(self) -> List[QuantumCircuit]:
                         (self.physical_qubits[1], self.physical_qubits[0]),
                     ]
 
+                qargs_supported = self.backend.target.qargs
                 instructions = []  # (op_name, qargs) for each element where qargs means qubit tuple
                 for qargs in qargs_patterns:
-                    for op_name in self.backend.target.operation_names_for_qargs(qargs):
-                        instructions.append((op_name, qargs))
+                    if qargs in qargs_supported:
+                        for op_name in self.backend.target.operation_names_for_qargs(qargs):
+                            instructions.append((op_name, qargs))
 
                 common_calibrations = defaultdict(dict)
                 for op_name, qargs in instructions:
diff --git a/releasenotes/notes/plugable-rb-clifford-synthesis-0e66c62fa3088fba.yaml b/releasenotes/notes/plugable-rb-clifford-synthesis-0e66c62fa3088fba.yaml
new file mode 100644
index 0000000000..1e5b36013d
--- /dev/null
+++ b/releasenotes/notes/plugable-rb-clifford-synthesis-0e66c62fa3088fba.yaml
@@ -0,0 +1,20 @@
+---
+features:
+  - |
+    Added a new experiment option ``clifford_synthesis_method`` to RB experiments,
+    e.g. :class:`~.StandardRB` and :class:`~.InterleavedRB` so that users can
+    plug in a custom Clifford synthesis algorithm used for generating RB circuits.
+    Such a plugin should be implemented as a ``HighLevelSynthesisPlugin``
+    (see :class:`~.RBDefaultCliffordSynthesis` for example).
+upgrade:
+  - |
+    Updated :class:`~.InterleavedRB` so that it only accepts ``interleaved_element``
+    consisting only of instructions supported by the backend of interest.
+fixes:
+  - |
+    Fixed a bug in circuit generation for three or more qubit RB where
+    sampled Cliffords may be changed during their circuits synthesis
+    (in the worst case, the resulting circuits may use qubits not in
+    ``physical_qubits``). See issue
+    `#1279 <https://github.com/Qiskit-Extensions/qiskit-experiments/issues/1279>`_
+    for additional details.
\ No newline at end of file
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 24c7719517..df9b3f3a08 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,4 +1,3 @@
-qiskit>=0.45.0
 black~=22.0
 fixtures
 stestr
diff --git a/requirements.txt b/requirements.txt
index 20f645b1a8..54ea5ea51c 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,6 +1,6 @@
 numpy>=1.17
 scipy>=1.4
-qiskit>=0.44
+qiskit>=0.45
 qiskit-ibm-experiment>=0.3.4
 matplotlib>=3.4
 uncertainties
diff --git a/setup.py b/setup.py
index 830fab8f9f..65ec0272a9 100755
--- a/setup.py
+++ b/setup.py
@@ -69,4 +69,9 @@
         "Source Code": "https://github.com/Qiskit-Extensions/qiskit-experiments",
     },
     zip_safe=False,
+    entry_points={
+        "qiskit.synthesis": [
+            "clifford.rb_default = qiskit_experiments.library.randomized_benchmarking.clifford_synthesis:RBDefaultCliffordSynthesis",
+        ],
+    },
 )
diff --git a/test/library/randomized_benchmarking/test_clifford_utils.py b/test/library/randomized_benchmarking/test_clifford_utils.py
index ddc11c9b94..37e68bf1e7 100644
--- a/test/library/randomized_benchmarking/test_clifford_utils.py
+++ b/test/library/randomized_benchmarking/test_clifford_utils.py
@@ -16,7 +16,7 @@
 from test.base import QiskitExperimentsTestCase
 
 import numpy as np
-from ddt import ddt
+from ddt import ddt, data
 from numpy.random import default_rng
 
 from qiskit import QuantumCircuit
@@ -32,7 +32,7 @@
     SXGate,
     RZGate,
 )
-from qiskit.quantum_info import Operator, Clifford
+from qiskit.quantum_info import Operator, Clifford, random_clifford
 from qiskit_experiments.library.randomized_benchmarking.clifford_utils import (
     CliffordUtils,
     num_from_1q_circuit,
@@ -45,6 +45,7 @@
     _layer_indices_from_num,
     _CLIFFORD_LAYER,
     _CLIFFORD_INVERSE_2Q,
+    _synthesize_clifford,
 )
 
 
@@ -221,3 +222,36 @@ def test_clifford_inverse_table(self):
         for lhs, rhs in enumerate(_CLIFFORD_INVERSE_2Q):
             c = compose_2q(lhs, rhs)
             self.assertEqual(c, 0)
+
+    @data(1, 2, 3, 4)
+    def test_clifford_synthesis_linear_connectivity(self, num_qubits):
+        """Check if clifford synthesis with linear connectivity does not change Clifford"""
+        basis_gates = tuple(["rz", "h", "cz"])
+        coupling_tuple = (
+            None if num_qubits == 1 else tuple((i, i + 1) for i in range(num_qubits - 1))
+        )
+        for seed in range(10):
+            expected = random_clifford(num_qubits=num_qubits, seed=seed)
+            circuit = _synthesize_clifford(expected, basis_gates, coupling_tuple)
+            synthesized = Clifford(circuit)
+            self.assertEqual(expected, synthesized)
+
+    @data(3, 4, 6)
+    def test_clifford_synthesis_non_linear_connectivity(self, num_qubits):
+        """Check if clifford synthesis with non-linear connectivity does not change Clifford"""
+        basis_gates = tuple(["rz", "sx", "cx"])
+        # star
+        coupling_tuple = tuple((0, i) for i in range(1, num_qubits))
+        for seed in range(5):
+            expected = random_clifford(num_qubits=num_qubits, seed=seed)
+            circuit = _synthesize_clifford(expected, basis_gates, coupling_tuple)
+            synthesized = Clifford(circuit)
+            self.assertEqual(expected, synthesized)
+
+        # cycle
+        coupling_tuple = tuple((i, (i + 1) % num_qubits) for i in range(num_qubits))
+        for seed in range(5):
+            expected = random_clifford(num_qubits=num_qubits, seed=seed)
+            circuit = _synthesize_clifford(expected, basis_gates, coupling_tuple)
+            synthesized = Clifford(circuit)
+            self.assertEqual(expected, synthesized)
diff --git a/test/library/randomized_benchmarking/test_interleaved_rb.py b/test/library/randomized_benchmarking/test_interleaved_rb.py
index 6ed28d0f53..1f7b567a23 100644
--- a/test/library/randomized_benchmarking/test_interleaved_rb.py
+++ b/test/library/randomized_benchmarking/test_interleaved_rb.py
@@ -16,10 +16,12 @@
 from test.library.randomized_benchmarking.mixin import RBTestMixin
 from ddt import ddt, data, unpack
 
+from qiskit import pulse
 from qiskit.circuit import Delay, QuantumCircuit, Parameter, Gate
 from qiskit.circuit.library import SXGate, CXGate, TGate, CZGate
 from qiskit.exceptions import QiskitError
 from qiskit.providers.fake_provider import FakeManila, FakeManilaV2, FakeWashington
+from qiskit.transpiler import InstructionProperties
 from qiskit_aer import AerSimulator
 from qiskit_aer.noise import NoiseModel, depolarizing_error
 from qiskit_experiments.library import randomized_benchmarking as rb
@@ -266,25 +268,44 @@ def test_interleaved_circuit_is_decomposed(self):
             self.assertTrue(all(not inst.operation.name.startswith("Clifford") for inst in qc))
 
     def test_interleaving_cnot_gate_with_non_supported_direction(self):
-        """Test if cx(0, 1) can be interleaved for backend that support only cx(1, 0)."""
+        """Test if fails to interleave cx(1, 2) for backend that support only cx(2, 1)."""
         my_backend = FakeManilaV2()
-        del my_backend.target["cx"][(0, 1)]  # make support only cx(1, 0)
+        del my_backend.target["cx"][(1, 2)]  # make support only cx(2, 1)
 
         exp = rb.InterleavedRB(
             interleaved_element=CXGate(),
-            physical_qubits=(0, 1),
+            physical_qubits=(1, 2),
             lengths=[3],
             num_samples=4,
             backend=my_backend,
             seed=1234,
         )
-        transpiled = exp._transpiled_circuits()
-        for qc in transpiled:
-            self.assertTrue(qc.count_ops().get("cx", 0) > 0)
-            expected_qubits = (qc.qubits[1], qc.qubits[0])
-            for inst in qc:
-                if inst.operation.name == "cx":
-                    self.assertEqual(inst.qubits, expected_qubits)
+        with self.assertRaises(QiskitError):
+            exp.circuits()
+
+    def test_interleaving_three_qubit_gate_with_calibration(self):
+        """Test if circuits for 3Q InterleavedRB contain custom calibrations supplied via target."""
+        my_backend = FakeManilaV2()
+        with pulse.build(my_backend) as custom_3q_sched:  # meaningless schedule
+            pulse.play(pulse.GaussianSquare(1600, 0.2, 64, 1300), pulse.drive_channel(0))
+
+        physical_qubits = (2, 1, 3)
+        custom_3q_gate = self.ThreeQubitGate()
+        my_backend.target.add_instruction(
+            custom_3q_gate, {physical_qubits: InstructionProperties(calibration=custom_3q_sched)}
+        )
+
+        exp = rb.InterleavedRB(
+            interleaved_element=custom_3q_gate,
+            physical_qubits=physical_qubits,
+            lengths=[3],
+            num_samples=1,
+            backend=my_backend,
+            seed=1234,
+        )
+        circuits = exp._transpiled_circuits()
+        qubits = tuple(circuits[0].qubits[q] for q in physical_qubits)
+        self.assertTrue(circuits[0].has_calibration_for((custom_3q_gate, qubits, [])))
 
 
 class TestRunInterleavedRB(QiskitExperimentsTestCase, RBTestMixin):
diff --git a/test/library/randomized_benchmarking/test_standard_rb.py b/test/library/randomized_benchmarking/test_standard_rb.py
index 71710292d0..fde71cc5e3 100644
--- a/test/library/randomized_benchmarking/test_standard_rb.py
+++ b/test/library/randomized_benchmarking/test_standard_rb.py
@@ -294,7 +294,7 @@ def test_two_qubit(self):
         self.assertAlmostEqual(epc.value.n, epc_expected, delta=0.3 * epc_expected)
 
     def test_three_qubit(self):
-        """Test two qubit RB. Use default basis gates."""
+        """Test three qubit RB. Use default basis gates."""
         exp = rb.StandardRB(
             physical_qubits=(0, 1, 2),
             lengths=list(range(1, 30, 3)),