diff --git a/cirq-core/cirq/transformers/dynamical_decoupling.py b/cirq-core/cirq/transformers/dynamical_decoupling.py
index e3caa7633c3..0c3ffb912d5 100644
--- a/cirq-core/cirq/transformers/dynamical_decoupling.py
+++ b/cirq-core/cirq/transformers/dynamical_decoupling.py
@@ -15,54 +15,57 @@
 """Transformer pass that adds dynamical decoupling operations to a circuit."""
 
 from functools import reduce
-from typing import Dict, Optional, Sequence, Tuple, Union
+from typing import Dict, Optional, Tuple, Union
+from itertools import cycle
 
 from cirq.transformers import transformer_api
+from cirq.transformers.analytical_decompositions import single_qubit_decompositions
+from cirq.transformers.analytical_decompositions import unitary_to_pauli_string
 import cirq
 import numpy as np
 
 
-def _repeat_sequence(
-    base_sequence: Sequence['cirq.Gate'], num_idle_moments: int
-) -> Sequence['cirq.Gate']:
-    """Returns the longest possible dynamical decoupling sequence."""
-    repeat_times = num_idle_moments // len(base_sequence)
-    return list(base_sequence) * repeat_times
-
-
-def _get_dd_sequence_from_schema_name(schema: str) -> Sequence['cirq.Gate']:
+def _get_dd_sequence_from_schema_name(schema: str) -> Tuple['cirq.Gate', ...]:
     """Gets dynamical decoupling sequence from a schema name."""
-    dd_sequence: Sequence['cirq.Gate']
     match schema:
+        case 'DEFAULT':
+            return (cirq.X, cirq.Y, cirq.X, cirq.Y)
         case 'XX_PAIR':
-            dd_sequence = (cirq.X, cirq.X)
+            return (cirq.X, cirq.X)
         case 'X_XINV':
-            dd_sequence = (cirq.X, cirq.X**-1)
+            return (cirq.X, cirq.X**-1)
         case 'YY_PAIR':
-            dd_sequence = (cirq.Y, cirq.Y)
+            return (cirq.Y, cirq.Y)
         case 'Y_YINV':
-            dd_sequence = (cirq.Y, cirq.Y**-1)
+            return (cirq.Y, cirq.Y**-1)
         case _:
             raise ValueError('Invalid schema name.')
-    return dd_sequence
 
 
-def _validate_dd_sequence(dd_sequence: Sequence['cirq.Gate']) -> None:
+def _pauli_up_to_global_phase(gate: 'cirq.Gate') -> Union['cirq.Pauli', None]:
+    for pauli_gate in [cirq.X, cirq.Y, cirq.Z]:
+        if cirq.equal_up_to_global_phase(gate, pauli_gate):
+            return pauli_gate
+    return None
+
+
+def _validate_dd_sequence(dd_sequence: Tuple['cirq.Gate', ...]) -> None:
     """Validates a given dynamical decoupling sequence.
 
     Args:
         dd_sequence: Input dynamical sequence to be validated.
 
-    Returns:
-        A tuple containing:
-            - is_valid (bool): True if the dd sequence is valid, False otherwise.
-            - error_message (str): An error message if the dd sequence is invalid, else None.
-
     Raises:
         ValueError: If dd_sequence is not valid.
     """
     if len(dd_sequence) < 2:
         raise ValueError('Invalid dynamical decoupling sequence. Expect more than one gates.')
+    for gate in dd_sequence:
+        if _pauli_up_to_global_phase(gate) is None:
+            raise ValueError(
+                'Dynamical decoupling sequence should only contain gates that are essentially'
+                ' Pauli gates.'
+            )
     matrices = [cirq.unitary(gate) for gate in dd_sequence]
     product = reduce(np.matmul, matrices)
 
@@ -73,14 +76,134 @@ def _validate_dd_sequence(dd_sequence: Sequence['cirq.Gate']) -> None:
         )
 
 
-def _parse_dd_sequence(schema: Union[str, Sequence['cirq.Gate']]) -> Sequence['cirq.Gate']:
+def _parse_dd_sequence(schema: Union[str, Tuple['cirq.Gate', ...]]) -> Tuple['cirq.Gate', ...]:
     """Parses and returns dynamical decoupling sequence from schema."""
     if isinstance(schema, str):
-        dd_sequence = _get_dd_sequence_from_schema_name(schema)
+        return _get_dd_sequence_from_schema_name(schema)
     else:
         _validate_dd_sequence(schema)
-        dd_sequence = schema
-    return dd_sequence
+        return schema
+
+
+def _is_single_qubit_operation(operation: 'cirq.Operation') -> bool:
+    if len(operation.qubits) != 1:
+        return False
+    return True
+
+
+def _is_single_qubit_gate_moment(moment: 'cirq.Moment') -> bool:
+    for operation in moment:
+        if not _is_single_qubit_operation(operation):
+            return False
+    return True
+
+
+def _is_clifford_moment(moment: 'cirq.Moment') -> bool:
+    for op in moment.operations:
+        if op.gate is not None and isinstance(op.gate, cirq.MeasurementGate):
+            return False
+        if not cirq.has_stabilizer_effect(op):
+            return False
+    return True
+
+
+def _get_clifford_pieces(circuit: 'cirq.AbstractCircuit') -> list[Tuple[int, int]]:
+    clifford_pieces: list[Tuple[int, int]] = []
+    left = 0
+    for moment_id, moment in enumerate(circuit):
+        if not _is_clifford_moment(moment):
+            clifford_pieces.append((left, moment_id))
+            left = moment_id + 1
+    if left < len(circuit):
+        clifford_pieces.append((left, len(circuit)))
+    return clifford_pieces
+
+
+def _is_insertable_moment(moment: 'cirq.Moment', single_qubit_gate_moments_only: bool) -> bool:
+    return _is_single_qubit_gate_moment(moment) or not single_qubit_gate_moments_only
+
+
+def _calc_pulled_through(
+    moment: 'cirq.Moment', input_pauli_ops: 'cirq.PauliString'
+) -> 'cirq.PauliString':
+    """Calculates the pulled_through after pulling through moment with the input.
+
+    We assume that the moment is Clifford here. Then, pulling through is essentially
+      decomposing a matrix into Pauli operations on each qubit.
+    """
+    pulled_through: 'cirq.PauliString' = cirq.PauliString()
+    for affected_q, combined_op_in_pauli in input_pauli_ops.items():
+        op_at_moment = moment.operation_at(affected_q)
+        if op_at_moment is None:
+            pulled_through *= combined_op_in_pauli.on(affected_q)
+            continue
+        prev_circuit = cirq.Circuit(cirq.Moment(op_at_moment))
+        new_circuit = cirq.Circuit(
+            cirq.Moment(combined_op_in_pauli.on(affected_q)), cirq.Moment(op_at_moment)
+        )
+        qubit_order = op_at_moment.qubits
+        pulled_through_pauli_ops = unitary_to_pauli_string(
+            prev_circuit.unitary(qubit_order=qubit_order)
+            @ new_circuit.unitary(qubit_order=qubit_order).conj().T
+        )
+        if pulled_through_pauli_ops is not None:
+            for qid, gate in enumerate(pulled_through_pauli_ops):
+                pulled_through *= gate.on(qubit_order[qid])
+    return pulled_through
+
+
+def _merge_pulled_through(
+    mutable_circuit: 'cirq.Circuit',
+    pulled_through: 'cirq.PauliString',
+    clifford_piece_range: Tuple[int, int],
+    single_qubit_gate_moments_only: bool,
+) -> 'cirq.PauliString':
+    """Merges pulled through Pauli gates into the last single-qubit gate operation or the insert it
+      into the first idle moment if idle moments exist.
+    Args:
+        mutable_circuit: Mutable circuit to transform.
+        pulled_through: Pauli gates to be merged.
+        clifford_piece_range: Specifies the [l, r) moments within which pulled-through gate merging
+          is to be performed.
+        single_qubit_gate_moments_only: If set True, dynamical decoupling operation will only be
+            added in single-qubit gate moments.
+
+    Returns:
+        The remaining pulled through operations after merging.
+    """
+    insert_intos: list[Tuple[int, 'cirq.Operation']] = []
+    batch_replaces: list[Tuple[int, 'cirq.Operation', 'cirq.Operation']] = []
+    remaining_pulled_through = pulled_through
+    for affected_q, combined_op_in_pauli in pulled_through.items():
+        moment_id = mutable_circuit.prev_moment_operating_on([affected_q], clifford_piece_range[1])
+        if moment_id is not None:
+            op = mutable_circuit.operation_at(affected_q, moment_id)
+            # Try to merge op into an existing single-qubit gate operation.
+            if op is not None and _is_single_qubit_operation(op):
+                updated_gate_mat = cirq.unitary(combined_op_in_pauli) @ cirq.unitary(op)
+                updated_gate: Optional['cirq.Gate'] = (
+                    single_qubit_decompositions.single_qubit_matrix_to_phxz(updated_gate_mat)
+                )
+                if updated_gate is None:
+                    # updated_gate is close to Identity.
+                    updated_gate = cirq.I
+                batch_replaces.append((moment_id, op, updated_gate.on(affected_q)))
+                remaining_pulled_through *= combined_op_in_pauli.on(affected_q)
+                continue
+            # Insert into the first empty moment for the qubit if such moment exists.
+            while moment_id < clifford_piece_range[1]:
+                if affected_q not in mutable_circuit.moments[
+                    moment_id
+                ].qubits and _is_insertable_moment(
+                    mutable_circuit.moments[moment_id], single_qubit_gate_moments_only
+                ):
+                    insert_intos.append((moment_id, combined_op_in_pauli.on(affected_q)))
+                    remaining_pulled_through *= combined_op_in_pauli.on(affected_q)
+                    break
+                moment_id += 1
+    mutable_circuit.batch_insert_into(insert_intos)
+    mutable_circuit.batch_replace(batch_replaces)
+    return remaining_pulled_through
 
 
 @transformer_api.transformer
@@ -88,10 +211,12 @@ def add_dynamical_decoupling(
     circuit: 'cirq.AbstractCircuit',
     *,
     context: Optional['cirq.TransformerContext'] = None,
-    schema: Union[str, Sequence['cirq.Gate']] = 'X_XINV',
+    schema: Union[str, Tuple['cirq.Gate', ...]] = 'DEFAULT',
+    single_qubit_gate_moments_only: bool = True,
 ) -> 'cirq.Circuit':
-    """Adds dynamical decoupling gate operations to idle moments of a given circuit.
-    This transformer preserves the moment structure of the circuit.
+    """Adds dynamical decoupling gate operations to a given circuit.
+    This transformer might add a new moment after each piece of Clifford moments, so the original
+      moment structure could change.
 
     Args:
           circuit: Input circuit to transform.
@@ -99,24 +224,70 @@ def add_dynamical_decoupling(
           schema: Dynamical decoupling schema name or a dynamical decoupling sequence.
             If a schema is specified, provided dynamical decouping sequence will be used.
             Otherwise, customized dynamical decoupling sequence will be applied.
+          single_qubit_gate_moments_only: If set True, dynamical decoupling operation will only be
+            added in single-qubit gate moments.
 
     Returns:
           A copy of the input circuit with dynamical decoupling operations.
     """
-    last_busy_moment_by_qubits: Dict['cirq.Qid', int] = {q: 0 for q in circuit.all_qubits()}
-    insert_into: list[Tuple[int, 'cirq.OP_TREE']] = []
+    base_dd_sequence: Tuple['cirq.Gate', ...] = _parse_dd_sequence(schema)
+    mutable_circuit = circuit.unfreeze(copy=True)
 
-    base_dd_sequence = _parse_dd_sequence(schema)
+    pauli_map: Dict['cirq.Gate', 'cirq.Pauli'] = {}
+    for gate in base_dd_sequence:
+        pauli_gate = _pauli_up_to_global_phase(gate)
+        if pauli_gate is not None:
+            pauli_map[gate] = pauli_gate
 
+    busy_moment_range_by_qubit: Dict['cirq.Qid', list[int]] = {
+        q: [len(circuit), -1] for q in circuit.all_qubits()
+    }
     for moment_id, moment in enumerate(circuit):
         for q in moment.qubits:
-            insert_gates = _repeat_sequence(
-                base_dd_sequence, num_idle_moments=moment_id - last_busy_moment_by_qubits[q] - 1
-            )
-            for idx, gate in enumerate(insert_gates):
-                insert_into.append((last_busy_moment_by_qubits[q] + idx + 1, gate.on(q)))
-            last_busy_moment_by_qubits[q] = moment_id
+            busy_moment_range_by_qubit[q][0] = min(busy_moment_range_by_qubit[q][0], moment_id)
+            busy_moment_range_by_qubit[q][1] = max(busy_moment_range_by_qubit[q][1], moment_id)
+    clifford_pieces = _get_clifford_pieces(circuit)
+
+    insert_intos: list[Tuple[int, 'cirq.Operation']] = []
+    insert_moments: list[Tuple[int, 'cirq.Moment']] = []
+    for l, r in clifford_pieces:  # [l, r)
+        # A PauliString stores the result of 'pulling' Pauli gates past each operations
+        # right before the current moment.
+        pulled_through: 'cirq.PauliString' = cirq.PauliString()
+        iter_by_qubits = {q: cycle(base_dd_sequence) for q in circuit.all_qubits()}
+
+        # Iterate over the Clifford piece.
+        for moment_id in range(l, r):
+            moment = circuit.moments[moment_id]
+
+            # Insert
+            if _is_insertable_moment(moment, single_qubit_gate_moments_only):
+                for q in circuit.all_qubits() - moment.qubits:
+                    if (
+                        busy_moment_range_by_qubit[q][0]
+                        < moment_id
+                        < busy_moment_range_by_qubit[q][1]
+                    ):
+                        insert_gate = next(iter_by_qubits[q])
+                        insert_intos.append((moment_id, insert_gate.on(q)))
+                        pulled_through *= pauli_map[insert_gate].on(q)
+
+            # Pull through
+            pulled_through = _calc_pulled_through(moment, pulled_through)
+
+        mutable_circuit.batch_insert_into(insert_intos)
+        insert_intos.clear()
+
+        pulled_through = _merge_pulled_through(
+            mutable_circuit, pulled_through, (l, r), single_qubit_gate_moments_only
+        )
+
+        # Insert a new moment if there are remaining pulled through operations.
+        new_moment_ops = []
+        for affected_q, combined_op_in_pauli in pulled_through.items():
+            new_moment_ops.append(combined_op_in_pauli.on(affected_q))
+        if len(new_moment_ops) != 0:
+            insert_moments.append((r, cirq.Moment(new_moment_ops)))
 
-    updated_circuit = circuit.unfreeze(copy=True)
-    updated_circuit.batch_insert_into(insert_into)
-    return updated_circuit
+    mutable_circuit.batch_insert(insert_moments)
+    return mutable_circuit
diff --git a/cirq-core/cirq/transformers/dynamical_decoupling_test.py b/cirq-core/cirq/transformers/dynamical_decoupling_test.py
index b7a17e28425..79a2bdd5317 100644
--- a/cirq-core/cirq/transformers/dynamical_decoupling_test.py
+++ b/cirq-core/cirq/transformers/dynamical_decoupling_test.py
@@ -12,26 +12,54 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-from typing import Sequence, Union
+from typing import Sequence, Tuple, Union
 import cirq
 from cirq import add_dynamical_decoupling
 import pytest
+import numpy as np
+
+
+def assert_sim_eq(circuit1: 'cirq.AbstractCircuit', circuit2: 'cirq.AbstractCircuit'):
+    # Simulate 2 circuits and compare final states.
+    sampler = cirq.Simulator(dtype=np.complex128)
+    psi0 = sampler.simulate(cirq.drop_terminal_measurements(circuit1)).final_state_vector
+    psi1 = sampler.simulate(cirq.drop_terminal_measurements(circuit2)).final_state_vector
+
+    assert np.isclose(np.abs(np.vdot(psi0, psi1)) ** 2, 1.0)
 
 
 def assert_dd(
-    input_circuit: cirq.Circuit,
-    expected_circuit: cirq.Circuit,
-    schema: Union[str, Sequence['cirq.Gate']],
+    input_circuit: 'cirq.AbstractCircuit',
+    expected_circuit: 'cirq.AbstractCircuit',
+    schema: Union[str, Tuple['cirq.Gate', ...]] = 'DEFAULT',
+    single_qubit_gate_moments_only: bool = True,
 ):
-    updated_circuit = add_dynamical_decoupling(input_circuit, schema=schema)
-    cirq.testing.assert_same_circuits(updated_circuit, expected_circuit)
+    transformed_circuit = add_dynamical_decoupling(
+        input_circuit, schema=schema, single_qubit_gate_moments_only=single_qubit_gate_moments_only
+    ).freeze()
+    cirq.testing.assert_same_circuits(transformed_circuit, expected_circuit)
+    cirq.testing.assert_circuits_have_same_unitary_given_final_permutation(
+        cirq.drop_terminal_measurements(input_circuit),
+        cirq.drop_terminal_measurements(transformed_circuit),
+        {q: q for q in input_circuit.all_qubits()},
+    )
+    assert_sim_eq(input_circuit, transformed_circuit)
 
 
-def test_no_insert_due_to_no_consecutive_moments():
+def test_no_insertion():
+    """Test case diagrams.
+    Input:
+    a: ───H───@───────
+              │
+    b: ───────X───H───
+    Output:
+    a: ───H───@───────
+              │
+    b: ───────X───H───
+    """
     a = cirq.NamedQubit('a')
     b = cirq.NamedQubit('b')
 
-    # No insertion as there is no room for a dd sequence.
     assert_dd(
         input_circuit=cirq.Circuit(
             cirq.Moment(cirq.H(a)), cirq.Moment(cirq.CNOT(a, b)), cirq.Moment(cirq.H(b))
@@ -40,6 +68,7 @@ def test_no_insert_due_to_no_consecutive_moments():
             cirq.Moment(cirq.H(a)), cirq.Moment(cirq.CNOT(a, b)), cirq.Moment(cirq.H(b))
         ),
         schema='XX_PAIR',
+        single_qubit_gate_moments_only=False,
     )
 
 
@@ -53,6 +82,14 @@ def test_no_insert_due_to_no_consecutive_moments():
     ],
 )
 def test_insert_provided_schema(schema: str, inserted_gates: Sequence['cirq.Gate']):
+    """Test case diagrams.
+    Input:
+    a: ───H───@───────────M───
+              │
+    b: ───────X───@───@───M───
+                  │   │
+    c: ───────────X───X───M───
+    """
     a = cirq.NamedQubit('a')
     b = cirq.NamedQubit('b')
     c = cirq.NamedQubit('c')
@@ -62,21 +99,35 @@ def test_insert_provided_schema(schema: str, inserted_gates: Sequence['cirq.Gate
         cirq.Moment(cirq.CNOT(a, b)),
         cirq.Moment(cirq.CNOT(b, c)),
         cirq.Moment(cirq.CNOT(b, c)),
-        cirq.Moment(cirq.measure_each(a, b, c)),
+        cirq.Moment([cirq.M(qubit) for qubit in [a, b, c]]),
     )
     expected_circuit = cirq.Circuit(
         cirq.Moment(cirq.H(a)),
         cirq.Moment(cirq.CNOT(a, b)),
         cirq.Moment(cirq.CNOT(b, c), inserted_gates[0](a)),
         cirq.Moment(cirq.CNOT(b, c), inserted_gates[1](a)),
-        cirq.Moment(cirq.measure_each(a, b, c)),
+        cirq.Moment([cirq.M(qubit) for qubit in [a, b, c]]),
     )
 
     # Insert one dynamical decoupling sequence in idle moments.
-    assert_dd(input_circuit, expected_circuit, schema=schema)
+    assert_dd(input_circuit, expected_circuit, schema=schema, single_qubit_gate_moments_only=False)
 
 
 def test_insert_by_customized_dd_sequence():
+    """Test case diagrams.
+        Input:
+    a: ───H───@───────────────────H───
+              │
+    b: ───────X───@───@───@───@───H───
+                  │   │   │   │
+    c: ───────────X───X───X───X───H───
+    Output:
+    a: ───H───@───X───X───Y───Y───H───
+              │
+    b: ───────X───@───@───@───@───H───
+                  │   │   │   │
+    c: ───────────X───X───X───X───H───
+    """
     a = cirq.NamedQubit('a')
     b = cirq.NamedQubit('b')
     c = cirq.NamedQubit('c')
@@ -89,7 +140,7 @@ def test_insert_by_customized_dd_sequence():
             cirq.Moment(cirq.CNOT(b, c)),
             cirq.Moment(cirq.CNOT(b, c)),
             cirq.Moment(cirq.CNOT(b, c)),
-            cirq.Moment(cirq.measure_each(a, b, c)),
+            cirq.Moment([cirq.H(qubit) for qubit in [a, b, c]]),
         ),
         expected_circuit=cirq.Circuit(
             cirq.Moment(cirq.H(a)),
@@ -98,9 +149,84 @@ def test_insert_by_customized_dd_sequence():
             cirq.Moment(cirq.CNOT(b, c), cirq.X(a)),
             cirq.Moment(cirq.CNOT(b, c), cirq.Y(a)),
             cirq.Moment(cirq.CNOT(b, c), cirq.Y(a)),
-            cirq.Moment(cirq.measure_each(a, b, c)),
+            cirq.Moment([cirq.H(qubit) for qubit in [a, b, c]]),
         ),
         schema=[cirq.X, cirq.X, cirq.Y, cirq.Y],
+        single_qubit_gate_moments_only=False,
+    )
+
+
+@pytest.mark.parametrize('single_qubit_gate_moments_only', [True, False])
+def test_pull_through_h_gate_case1(single_qubit_gate_moments_only: bool):
+    """Test case diagrams.
+    Input:
+    a: ───H───────H───────@───
+                          │
+    b: ───H───H───H───H───X───
+    Output:
+    a: ───H───X───H───X───@───Y───
+                          │
+    b: ───H───H───H───H───X───X───
+    """
+    a = cirq.NamedQubit('a')
+    b = cirq.NamedQubit('b')
+
+    assert_dd(
+        input_circuit=cirq.Circuit(
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b)),
+            cirq.Moment(cirq.CNOT(a, b)),
+        ),
+        expected_circuit=cirq.Circuit(
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b), cirq.X(a)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b), cirq.X(a)),
+            cirq.Moment(cirq.CNOT(a, b)),
+            cirq.Moment(cirq.Y(a), cirq.X(b)),
+        ),
+        schema="XX_PAIR",
+        single_qubit_gate_moments_only=single_qubit_gate_moments_only,
+    )
+
+
+@pytest.mark.parametrize('single_qubit_gate_moments_only', [True, False])
+def test_pull_through_h_gate_case2(single_qubit_gate_moments_only: bool):
+    """Test case diagrams.
+    Input:
+    a: ───H───────H───────H───
+
+    b: ───H───H───H───H───H───
+    Output:
+    a: ───H───X───H───X───PhXZ(a=0.5,x=0.5,z=1)───
+
+    b: ───H───H───H───H───H───────────────────────
+    """
+    a = cirq.NamedQubit('a')
+    b = cirq.NamedQubit('b')
+
+    assert_dd(
+        input_circuit=cirq.Circuit(
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+        ),
+        expected_circuit=cirq.Circuit(
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b), cirq.X(a)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b), cirq.X(a)),
+            cirq.Moment(
+                cirq.PhasedXZGate(axis_phase_exponent=0.5, x_exponent=0.5, z_exponent=1).on(a),
+                cirq.H(b),
+            ),
+        ),
+        schema="XX_PAIR",
+        single_qubit_gate_moments_only=single_qubit_gate_moments_only,
     )
 
 
@@ -110,14 +236,393 @@ def test_insert_by_customized_dd_sequence():
         ('INVALID_SCHEMA', 'Invalid schema name.'),
         ([cirq.X], 'Invalid dynamical decoupling sequence. Expect more than one gates.'),
         (
-            [cirq.X, cirq.H],
+            [cirq.X, cirq.Y],
             'Invalid dynamical decoupling sequence. Expect sequence production equals identity'
             ' up to a global phase, got',
         ),
+        (
+            [cirq.H, cirq.H],
+            'Dynamical decoupling sequence should only contain gates that are essentially'
+            ' Pauli gates.',
+        ),
     ],
 )
-def test_invalid_dd_schema(schema: Union[str, Sequence['cirq.Gate']], error_msg_regex):
+def test_invalid_dd_schema(schema: Union[str, Tuple['cirq.Gate', ...]], error_msg_regex):
     a = cirq.NamedQubit('a')
     input_circuit = cirq.Circuit(cirq.H(a))
     with pytest.raises(ValueError, match=error_msg_regex):
-        add_dynamical_decoupling(input_circuit, schema=schema)
+        add_dynamical_decoupling(input_circuit, schema=schema, single_qubit_gate_moments_only=False)
+
+
+def test_single_qubit_gate_moments_only_no_updates_succeeds():
+    qubits = cirq.LineQubit.range(9)
+    input_circuit = cirq.Circuit(
+        cirq.Moment([cirq.H(qubits[i]) for i in [3, 4, 5]]),
+        cirq.Moment(cirq.CZ(*qubits[4:6])),
+        cirq.Moment(cirq.CZ(*qubits[3:5])),
+        cirq.Moment([cirq.H(qubits[i]) for i in [2, 3, 5, 6]]),
+        cirq.Moment(cirq.CZ(*qubits[2:4]), cirq.CNOT(*qubits[5:7])),
+        cirq.Moment([cirq.H(qubits[i]) for i in [1, 2, 6, 7]]),
+        cirq.Moment(cirq.CZ(*qubits[1:3]), cirq.CNOT(*qubits[6:8])),
+        cirq.Moment([cirq.H(qubits[i]) for i in [0, 1, 7, 8]]),
+        cirq.Moment(cirq.CZ(*qubits[0:2]), cirq.CNOT(*qubits[7:])),
+    )
+    add_dynamical_decoupling(input_circuit, schema='X_XINV', single_qubit_gate_moments_only=True)
+
+
+def test_scattered_circuit():
+    """Test case diagrams.
+    Input:
+    0: ───────────────────────────────H───@───H───
+                                          │
+    1: ───────────────────────H───@───H───@───H───
+                                  │
+    2: ───────────────H───@───H───@───────────H───
+                          │
+    3: ───H───────@───H───@───────────────────H───
+                  │
+    4: ───H───@───@───────────────────────────H───
+              │
+    5: ───H───@───────H───@───────────────────H───
+                          │
+    6: ───────────────H───@───H───@───────────H───
+                                  │
+    7: ───────────────────────H───@───H───@───H───
+                                          │
+    8: ───────────────────────────────H───@───H───
+
+    Output (single_qubit_gate_moment_only_on):
+    0: ───────────────────────────────H───@───H────────────────────────
+                                          │
+    1: ───────────────────────H───@───H───@───H────────────────────────
+                                  │
+    2: ───────────────H───@───H───@───X───────PhXZ(a=-0.5,x=0.5,z=0)───
+                          │
+    3: ───H───────@───H───@───X───────Y───────PhXZ(a=0.5,x=0.5,z=0)────
+                  │
+    4: ───H───@───@───X───────Y───────X───────PhXZ(a=0.5,x=0.5,z=-1)───
+              │
+    5: ───H───@───────H───@───X───────Y───────PhXZ(a=0.5,x=0.5,z=0)────
+                          │
+    6: ───────────────H───@───H───@───X───────PhXZ(a=-0.5,x=0.5,z=0)───
+                                  │
+    7: ───────────────────────H───@───H───@───H────────────────────────
+                                          │
+    8: ───────────────────────────────H───@───H────────────────────────
+
+    Output (single_qubit_gate_moment_only_off):
+    0: ───────────────────────────────H───@───H───────────────────────
+                                          │
+    1: ───────────────────────H───@───H───@───H───────────────────────
+                                  │
+    2: ───────────────H───@───H───@───X───Y───PhXZ(a=0.5,x=0.5,z=0)───
+                          │
+    3: ───H───X───@───H───@───Y───X───Y───X───PhXZ(a=0.5,x=0.5,z=0)───
+                  │
+    4: ───H───@───@───X───Y───X───Y───X───Y───H───────────────────────
+              │
+    5: ───H───@───X───H───@───Y───X───Y───X───PhXZ(a=0.5,x=0.5,z=0)───
+                          │
+    6: ───────────────H───@───H───@───X───Y───PhXZ(a=0.5,x=0.5,z=0)───
+                                  │
+    7: ───────────────────────H───@───H───@───H───────────────────────
+                                          │
+    8: ───────────────────────────────H───@───H───────────────────────
+    """
+    qubits = cirq.LineQubit.range(9)
+    input_circuit = cirq.Circuit(
+        cirq.Moment([cirq.H(qubits[i]) for i in [3, 4, 5]]),
+        cirq.Moment(cirq.CZ(*qubits[4:6])),
+        cirq.Moment(cirq.CZ(*qubits[3:5])),
+        cirq.Moment([cirq.H(qubits[i]) for i in [2, 3, 5, 6]]),
+        cirq.Moment(cirq.CZ(*qubits[2:4]), cirq.CZ(*qubits[5:7])),
+        cirq.Moment([cirq.H(qubits[i]) for i in [1, 2, 6, 7]]),
+        cirq.Moment(cirq.CZ(*qubits[1:3]), cirq.CZ(*qubits[6:8])),
+        cirq.Moment([cirq.H(qubits[i]) for i in [0, 1, 7, 8]]),
+        cirq.Moment(cirq.CZ(*qubits[0:2]), cirq.CZ(*qubits[7:])),
+        cirq.Moment([cirq.H(q) for q in qubits]),
+    )
+    expected_circuit_single_qubit_gate_on = cirq.Circuit(
+        cirq.Moment([cirq.H(qubits[i]) for i in [3, 4, 5]]),
+        cirq.Moment(cirq.CZ(*qubits[4:6])),
+        cirq.Moment(cirq.CZ(*qubits[3:5])),
+        cirq.Moment([cirq.H(qubits[i]) for i in [2, 3, 5, 6]] + [cirq.X(qubits[4])]),
+        cirq.Moment(cirq.CZ(*qubits[2:4]), cirq.CZ(*qubits[5:7])),
+        cirq.Moment(
+            [cirq.H(qubits[i]) for i in [1, 2, 6, 7]]
+            + [cirq.X(qubits[i]) for i in [3, 5]]
+            + [cirq.Y(qubits[4])]
+        ),
+        cirq.Moment(cirq.CZ(*qubits[1:3]), cirq.CZ(*qubits[6:8])),
+        cirq.Moment(
+            [cirq.H(qubits[i]) for i in [0, 1, 7, 8]]
+            + [cirq.X(qubits[i]) for i in [2, 4, 6]]
+            + [cirq.Y(qubits[i]) for i in [3, 5]]
+        ),
+        cirq.Moment(cirq.CZ(*qubits[0:2]), cirq.CZ(*qubits[7:])),
+        cirq.Moment(
+            [cirq.H(qubits[i]) for i in [0, 1, 7, 8]]
+            + [
+                cirq.PhasedXZGate(axis_phase_exponent=-0.5, x_exponent=0.5, z_exponent=0).on(
+                    qubits[2]
+                ),
+                cirq.PhasedXZGate(axis_phase_exponent=0.5, x_exponent=0.5, z_exponent=0).on(
+                    qubits[3]
+                ),
+                cirq.PhasedXZGate(axis_phase_exponent=0.5, x_exponent=0.5, z_exponent=-1).on(
+                    qubits[4]
+                ),
+                cirq.PhasedXZGate(axis_phase_exponent=0.5, x_exponent=0.5, z_exponent=0).on(
+                    qubits[5]
+                ),
+                cirq.PhasedXZGate(axis_phase_exponent=-0.5, x_exponent=0.5, z_exponent=0).on(
+                    qubits[6]
+                ),
+            ]
+        ),
+    )
+    expected_circuit_single_qubit_gates_off = cirq.Circuit(
+        cirq.Moment([cirq.H(qubits[i]) for i in [3, 4, 5]]),
+        cirq.Moment(cirq.CZ(*qubits[4:6]), cirq.X(qubits[3])),
+        cirq.Moment(cirq.CZ(*qubits[3:5]), cirq.X(qubits[5])),
+        cirq.Moment([cirq.H(qubits[i]) for i in [2, 3, 5, 6]] + [cirq.X(qubits[i]) for i in [4]]),
+        cirq.Moment(cirq.CZ(*qubits[2:4]), cirq.CZ(*qubits[5:7]), cirq.Y(qubits[4])),
+        cirq.Moment(
+            [cirq.H(qubits[i]) for i in [1, 2, 6, 7]]
+            + [cirq.Y(qubits[i]) for i in [3, 5]]
+            + [cirq.X(qubits[4])]
+        ),
+        cirq.Moment(
+            [cirq.CZ(*qubits[1:3]), cirq.CZ(*qubits[6:8])]
+            + [cirq.X(qubits[i]) for i in [3, 5]]
+            + [cirq.Y(qubits[4])]
+        ),
+        cirq.Moment(
+            [cirq.H(qubits[i]) for i in [0, 1, 7, 8]]
+            + [cirq.X(qubits[i]) for i in [2, 4, 6]]
+            + [cirq.Y(qubits[i]) for i in [3, 5]]
+        ),
+        cirq.Moment(
+            [cirq.CZ(*qubits[0:2]), cirq.CZ(*qubits[7:])]
+            + [cirq.X(qubits[i]) for i in [3, 5]]
+            + [cirq.Y(qubits[i]) for i in [2, 4, 6]]
+        ),
+        cirq.Moment(
+            [cirq.H(qubits[i]) for i in [0, 1, 4, 7, 8]]
+            + [
+                cirq.PhasedXZGate(axis_phase_exponent=0.5, x_exponent=0.5, z_exponent=0).on(
+                    qubits[i]
+                )
+                for i in [2, 3, 5, 6]
+            ]
+        ),
+    )
+    assert_dd(
+        input_circuit,
+        expected_circuit_single_qubit_gate_on,
+        schema='DEFAULT',
+        single_qubit_gate_moments_only=True,
+    )
+    assert_dd(
+        input_circuit,
+        expected_circuit_single_qubit_gates_off,
+        schema='DEFAULT',
+        single_qubit_gate_moments_only=False,
+    )
+
+
+def test_scattered_circuit2():
+    """Test case diagrams.
+    Input:
+    0: ───────────────────@───
+                          │
+    1: ───────────────@───@───
+                      │
+    2: ───────────@───@───────
+                  │
+    3: ───────@───@───────────
+              │
+    4: ───@───@───────────────
+          │
+    5: ───@───────@───────────
+                  │
+    6: ───────────@───@───────
+                      │
+    7: ───────────────@───@───
+                          │
+    8: ───────────────────@───
+    Output:
+    0: ───────────────────@───
+                          │
+    1: ───────────────@───@───
+                      │
+    2: ───────────@───@───────
+                  │
+    3: ───────@───@───────────
+              │
+    4: ───@───@───────────────
+          │
+    5: ───@───X───@───X───────
+                  │
+    6: ───────────@───@───Z───
+                      │
+    7: ───────────────@───@───
+                          │
+    8: ───────────────────@───
+    """
+    qubits = cirq.LineQubit.range(9)
+    assert_dd(
+        input_circuit=cirq.Circuit(
+            cirq.Moment(cirq.CZ(*qubits[4:6])),
+            cirq.Moment(cirq.CZ(*qubits[3:5])),
+            cirq.Moment(cirq.CZ(*qubits[2:4]), cirq.CZ(*qubits[5:7])),
+            cirq.Moment(cirq.CZ(*qubits[1:3]), cirq.CZ(*qubits[6:8])),
+            cirq.Moment(cirq.CZ(*qubits[0:2]), cirq.CZ(*qubits[7:])),
+        ),
+        expected_circuit=cirq.Circuit(
+            cirq.Moment(cirq.CZ(*qubits[4:6])),
+            cirq.Moment(cirq.CZ(*qubits[3:5]), cirq.X(qubits[5])),
+            cirq.Moment(cirq.CZ(*qubits[2:4]), cirq.CZ(*qubits[5:7])),
+            cirq.Moment(cirq.CZ(*qubits[1:3]), cirq.CZ(*qubits[6:8]), cirq.X(qubits[5])),
+            cirq.Moment(cirq.CZ(*qubits[0:2]), cirq.CZ(*qubits[7:]), cirq.Z(qubits[6])),
+        ),
+        schema="XX_PAIR",
+        single_qubit_gate_moments_only=False,
+    )
+
+
+def test_pull_through_chain():
+    """Test case diagrams.
+    Input:
+    0: ───X───────×───────────X───
+                  │
+    1: ───────Y───×───×───────X───
+                      │
+    2: ───────────────×───×───X───
+                          │
+    3: ───────────────────×───X───
+    Output:
+    0: ───X───X───×───X───X───X───
+                  │
+    1: ───────Y───×───×───X───I───
+                      │
+    2: ───────────────×───×───X───
+                          │
+    3: ───────────────────×───I───
+    """
+    qubits = cirq.LineQubit.range(4)
+    assert_dd(
+        input_circuit=cirq.Circuit(
+            cirq.Moment(cirq.X(qubits[0])),
+            cirq.Moment(cirq.Y(qubits[1])),
+            cirq.Moment(cirq.SWAP(*qubits[0:2])),
+            cirq.Moment(cirq.SWAP(*qubits[1:3])),
+            cirq.Moment(cirq.SWAP(*qubits[2:4])),
+            cirq.Moment([cirq.X(qubits[i]) for i in range(4)]),
+        ),
+        expected_circuit=cirq.Circuit(
+            cirq.Moment(cirq.X(qubits[0])),
+            cirq.Moment(cirq.Y(qubits[1]), cirq.X(qubits[0])),
+            cirq.Moment(cirq.SWAP(*qubits[0:2])),
+            cirq.Moment([cirq.SWAP(*qubits[1:3])] + [cirq.X(qubits[0])]),
+            cirq.Moment([cirq.SWAP(*qubits[2:4])] + [cirq.X(qubits[0]), cirq.X(qubits[1])]),
+            cirq.Moment(cirq.X(qubits[0]), cirq.I(qubits[1]), cirq.X(qubits[2]), cirq.I(qubits[3])),
+        ),
+        schema='XX_PAIR',
+        single_qubit_gate_moments_only=False,
+    )
+
+
+def test_multiple_clifford_pieces():
+    """Test case diagrams.
+    Input:
+    a: ───H───────H───────@───────────H───────H───
+                          │
+    b: ───H───H───H───H───@^0.5───H───H───H───H───
+    Output:
+    a: ───H───X───H───PhXZ(a=0.5,x=0,z=-1)───@───────X───H───X───PhXZ(a=0.5,x=0.5,z=-1)───
+                                             │
+    b: ───H───H───H───H──────────────────────@^0.5───H───H───H───H────────────────────────
+    """
+    a = cirq.NamedQubit('a')
+    b = cirq.NamedQubit('b')
+    assert_dd(
+        input_circuit=cirq.Circuit(
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b)),
+            cirq.Moment(cirq.CZPowGate(exponent=0.5).on(a, b)),
+            cirq.Moment(cirq.H(b)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+        ),
+        expected_circuit=cirq.Circuit(
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b), cirq.X(a)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(
+                cirq.H(b),
+                cirq.PhasedXZGate(axis_phase_exponent=0.5, x_exponent=0, z_exponent=-1).on(a),
+            ),
+            cirq.Moment(cirq.CZPowGate(exponent=0.5).on(a, b)),
+            cirq.Moment(cirq.H(b), cirq.X(a)),
+            cirq.Moment(cirq.H(a), cirq.H(b)),
+            cirq.Moment(cirq.H(b), cirq.X(a)),
+            cirq.Moment(
+                cirq.H(b),
+                cirq.PhasedXZGate(axis_phase_exponent=0.5, x_exponent=0.5, z_exponent=-1).on(a),
+            ),
+        ),
+        schema="XX_PAIR",
+    )
+
+
+def test_with_non_clifford_measurements():
+    """Test case diagrams.
+    Input:
+    0: ───────────H───@───H───M───
+                      │
+    1: ───H───@───────@───────M───
+              │
+    2: ───H───@───H───@───────M───
+                      │
+    3: ───────────H───@───H───M───
+    Output:
+    0: ───────────H───@───PhXZ(a=0.5,x=0.5,z=0)───M───
+                      │
+    1: ───H───@───X───@───X───────────────────────M───
+              │
+    2: ───H───@───H───@───I───────────────────────M───
+                      │
+    3: ───────────H───@───H───────────────────────M───
+    """
+    qubits = cirq.LineQubit.range(4)
+    assert_dd(
+        input_circuit=cirq.Circuit(
+            cirq.Moment([cirq.H(qubits[i]) for i in [1, 2]]),
+            cirq.Moment(cirq.CZ(*qubits[1:3])),
+            cirq.Moment([cirq.H(qubits[i]) for i in [0, 2, 3]]),
+            cirq.Moment(cirq.CZ(*qubits[0:2]), cirq.CZ(*qubits[2:])),
+            cirq.Moment([cirq.H(qubits[i]) for i in [0, 3]]),
+            cirq.Moment([cirq.M(qubits[i]) for i in [0, 1, 2, 3]]),
+        ),
+        expected_circuit=cirq.Circuit(
+            cirq.Moment([cirq.H(qubits[i]) for i in [1, 2]]),
+            cirq.Moment(cirq.CZ(*qubits[1:3])),
+            cirq.Moment([cirq.H(qubits[i]) for i in [0, 2, 3]] + [cirq.X(qubits[1])]),
+            cirq.Moment(cirq.CZ(*qubits[0:2]), cirq.CZ(*qubits[2:])),
+            cirq.Moment(
+                cirq.H(qubits[3]),
+                cirq.I(qubits[2]),
+                cirq.X(qubits[1]),
+                cirq.PhasedXZGate(axis_phase_exponent=0.5, x_exponent=0.5, z_exponent=0).on(
+                    qubits[0]
+                ),
+            ),
+            cirq.Moment([cirq.M(qubits[i]) for i in [0, 1, 2, 3]]),
+        ),
+        schema="XX_PAIR",
+        single_qubit_gate_moments_only=True,
+    )