diff --git a/src/tsim/circuit.py b/src/tsim/circuit.py
index 09210567..422d32fc 100644
--- a/src/tsim/circuit.py
+++ b/src/tsim/circuit.py
@@ -183,6 +183,49 @@ def append(
                 tag = f"U3(theta={theta}*pi, phi={phi}*pi, lambda={lam}*pi)"
                 name = "I"
                 arg = None
+            elif name in ("R_XX", "R_YY", "R_ZZ"):
+                if arg is None:
+                    raise ValueError(f"For {name} gates, an angle must be provided.")
+                args = list(arg) if isinstance(arg, Iterable) else [arg]
+                if len(args) != 1:
+                    raise ValueError(
+                        f"For {name} gates, a single angle must be provided."
+                    )
+                t_list: list[int] = (
+                    [targets] if isinstance(targets, int) else list(targets)  # type: ignore[arg-type]
+                )
+                if len(t_list) != 2:
+                    raise ValueError(
+                        f"For {name} gates, exactly two target qubits must be provided."
+                    )
+                q0, q1 = (
+                    t.value if isinstance(t, stim.GateTarget) else t for t in t_list
+                )
+                if q0 == q1:
+                    raise ValueError(
+                        f"Duplicate target qubits in {name}: both targets are qubit {q0}."
+                    )
+                pauli = name[2]  # 'X', 'Y', or 'Z'
+                target_fn = {
+                    "X": stim.target_x,
+                    "Y": stim.target_y,
+                    "Z": stim.target_z,
+                }[pauli]
+                tag = f"{name}(theta={args[0]}*pi)"
+                targets = [target_fn(q0), stim.target_combiner(), target_fn(q1)]
+                name = "SPP"
+                arg = None
+            elif name == "R_PAULI":
+                if arg is None:
+                    raise ValueError("For R_PAULI gates, an angle must be provided.")
+                args = list(arg) if isinstance(arg, Iterable) else [arg]
+                if len(args) != 1:
+                    raise ValueError(
+                        "For R_PAULI gates, a single angle must be provided."
+                    )
+                tag = f"R_PAULI(theta={args[0]}*pi)"
+                name = "SPP"
+                arg = None
 
             self._stim_circ.append(name=name, targets=targets, arg=arg, tag=tag)  # type: ignore
         else:
@@ -806,6 +849,25 @@ def fix_tags(circuit: stim.Circuit) -> stim.Circuit:
                         result.append("I", targets, args, tag=new_tag)
                         continue
 
+                # Stim flips SPP↔SPP_DAG during inverse.  For parametric
+                # Pauli rotations the angle already encodes direction, so
+                # undo the flip and negate θ instead.
+                if name in ("SPP", "SPP_DAG") and tag and tag != "T":
+                    parsed = parse_parametric_tag(instr)
+                    if parsed is not None:
+                        gate_name, params = parsed
+                        if name == "SPP_DAG":
+                            # Stim flipped SPP → SPP_DAG; undo flip by negating θ.
+                            theta = float(-params["theta"])
+                            new_tag = f"{gate_name}(theta={theta}*pi)"
+                            result.append(
+                                "SPP", instr.targets_copy(), args, tag=new_tag
+                            )
+                        else:
+                            # Stim flipped SPP_DAG → SPP; tag already correct.
+                            result.append(instr)
+                        continue
+
                 result.append(instr)
             return result
 
diff --git a/src/tsim/core/instructions.py b/src/tsim/core/instructions.py
index 0e62d734..48f84b74 100644
--- a/src/tsim/core/instructions.py
+++ b/src/tsim/core/instructions.py
@@ -963,6 +963,31 @@ def _pauli_product_phase(
             s(b, qubit)
 
 
+def r_pauli(
+    b: GraphRepresentation,
+    paulis: list[tuple[Literal["X", "Y", "Z"], int]],
+    phase: Fraction,
+) -> None:
+    """Apply exp(-i * phase * pi/2 * P) for an arbitrary Pauli product P.
+
+    Generalises R_X / R_Y / R_Z to multi-qubit Pauli strings.
+    ``R_XX``, ``R_YY``, ``R_ZZ`` are two-qubit special cases.
+
+    Args:
+        b: The graph representation to modify.
+        paulis: List of (pauli_type, qubit) pairs defining the Pauli product P.
+        phase: Rotation angle in units of π  (α means exp(-i α π/2 P)).
+
+    """
+    _pauli_product_phase(
+        b,
+        paulis,
+        lambda b_, q: r_z(b_, q, phase),
+        lambda b_, q: r_z(b_, q, -phase),
+        dagger=False,
+    )
+
+
 def spp(
     b: GraphRepresentation,
     paulis: list[tuple[Literal["X", "Y", "Z"], int]],
diff --git a/src/tsim/core/parse.py b/src/tsim/core/parse.py
index 0e74c712..0b4c2afc 100644
--- a/src/tsim/core/parse.py
+++ b/src/tsim/core/parse.py
@@ -16,6 +16,7 @@
     mpad,
     mpp,
     observable_include,
+    r_pauli,
     r_x,
     r_y,
     r_z,
@@ -29,6 +30,10 @@
     "R_X": frozenset({"theta"}),
     "R_Y": frozenset({"theta"}),
     "R_Z": frozenset({"theta"}),
+    "R_XX": frozenset({"theta"}),
+    "R_YY": frozenset({"theta"}),
+    "R_ZZ": frozenset({"theta"}),
+    "R_PAULI": frozenset({"theta"}),
     "U3": frozenset({"theta", "phi", "lambda"}),
 }
 
@@ -38,7 +43,7 @@ def parse_parametric_tag(
 ) -> tuple[str, dict[str, Fraction]] | None:
     """Parse the parametric tag on an instruction (e.g. ``I[R_Z(theta=0.3*pi)]``).
 
-    Supports gates: R_Z, R_X, R_Y, U3.
+    Supports gates: R_Z, R_X, R_Y, R_XX, R_YY, R_ZZ, R_PAULI, U3.
 
     Args:
         instruction: The stim instruction whose tag will be parsed.
@@ -225,6 +230,17 @@ def parse_stim_circuit(
             for paulis, invert in _iter_pauli_products(instruction):
                 mpp(b, paulis, invert, p=p)
             continue
+        if name in ("SPP", "SPP_DAG") and instruction.tag and instruction.tag != "T":
+            result = parse_parametric_tag(instruction)
+            if result is not None:
+                gate_name, params = result
+                is_dag = name == "SPP_DAG"
+                for paulis, invert in _iter_pauli_products(instruction):
+                    phase = params["theta"]
+                    if is_dag ^ invert:
+                        phase = -phase
+                    r_pauli(b, paulis, phase)
+                continue
         if name in ("SPP", "SPP_DAG") and instruction.tag == "T":
             is_dag = name == "SPP_DAG"
             for paulis, invert in _iter_pauli_products(instruction):
diff --git a/src/tsim/utils/clifford.py b/src/tsim/utils/clifford.py
index f27cc06f..c5c82ad3 100644
--- a/src/tsim/utils/clifford.py
+++ b/src/tsim/utils/clifford.py
@@ -6,7 +6,7 @@
 
 import stim
 
-from tsim.core.parse import parse_parametric_tag
+from tsim.core.parse import _iter_pauli_products, parse_parametric_tag
 
 # Clifford decompositions for U3(θ, φ, λ) = R_Z(φ) · R_Y(θ) · R_Z(λ).
 # Keys: (θ_idx, φ_idx, λ_idx) where each index ∈ {0,1,2,3} is the angle in half-pi units.
@@ -84,6 +84,20 @@ def parametric_to_clifford_gates(
         table = {"R_Z": RZ_CLIFFORD, "R_X": RX_CLIFFORD, "R_Y": RY_CLIFFORD}[gate_name]
         return [table[idx]]
 
+    if gate_name in ("R_XX", "R_YY", "R_ZZ", "R_PAULI"):
+        idx = _to_half_pi_index(params["theta"])
+        if idx is None:
+            return None
+        # Clifford-angle Pauli rotations map to SPP / SPP_DAG / identity
+        # idx=0 → I, idx=1 → SPP, idx=2 → SPP·SPP, idx=3 → SPP_DAG
+        _SPP_CLIFFORD: dict[int, list[str]] = {
+            0: [],
+            1: ["SPP"],
+            2: ["SPP", "SPP"],
+            3: ["SPP_DAG"],
+        }
+        return _SPP_CLIFFORD[idx]
+
     if gate_name == "U3":
         theta_idx = _to_half_pi_index(params["theta"])
         phi_idx = _to_half_pi_index(params["phi"])
@@ -119,6 +133,13 @@ def is_half_pi_multiple(phase: Fraction) -> bool:
         if instr.name in ["S", "S_DAG", "SPP", "SPP_DAG"] and instr.tag == "T":
             return False
 
+        if instr.name in ["SPP", "SPP_DAG"] and instr.tag and instr.tag != "T":
+            result = parse_parametric_tag(instr)
+            if result is not None:
+                _, params = result
+                if not is_half_pi_multiple(params["theta"]):
+                    return False
+
         if instr.name == "I" and instr.tag:
             result = parse_parametric_tag(instr)
             if result is None:
@@ -154,6 +175,11 @@ def expand_clifford_rotations(source: stim.Circuit) -> stim.Circuit:
                 )
             )
             continue
+        spp_exp = _try_spp_clifford_expansion(instr)
+        if spp_exp is not None:
+            for gate_name, gate_targets in spp_exp:
+                out.append(gate_name, gate_targets, [])
+            continue
         expansion = _try_clifford_expansion(instr)
         if expansion is not None:
             gates, targets = expansion
@@ -164,6 +190,38 @@ def expand_clifford_rotations(source: stim.Circuit) -> stim.Circuit:
     return out
 
 
+def _try_spp_clifford_expansion(
+    instr: stim.CircuitInstruction,
+) -> list[tuple[str, list[stim.GateTarget]]] | None:
+    """Try to expand a tagged ``SPP`` instruction into Clifford SPP/SPP_DAG.
+
+    Returns:
+        List of ``(gate_name, targets)`` pairs, or ``None`` if the instruction
+        is not an expandable parametric Pauli rotation.
+
+    """
+    if instr.name not in ("SPP", "SPP_DAG") or not instr.tag or instr.tag == "T":
+        return None
+
+    parsed = parse_parametric_tag(instr)
+    if parsed is None:
+        return None
+
+    _, params = parsed
+    effective_params = dict(params)
+    is_dag = instr.name == "SPP_DAG"
+    for _, invert in _iter_pauli_products(instr):
+        if is_dag ^ invert:
+            effective_params["theta"] = -effective_params["theta"]
+        break
+    gates = parametric_to_clifford_gates(parsed[0], effective_params)
+    if gates is None:
+        return None
+
+    targets = instr.targets_copy()
+    return [(g, targets) for g in gates] if gates else []
+
+
 def _try_clifford_expansion(
     instr: stim.CircuitInstruction,
 ) -> tuple[list[str], list[int]] | None:
diff --git a/src/tsim/utils/program_text.py b/src/tsim/utils/program_text.py
index ed8f6f41..7239dbde 100644
--- a/src/tsim/utils/program_text.py
+++ b/src/tsim/utils/program_text.py
@@ -5,9 +5,11 @@
 # Matches valid numeric literals including scientific notation (e.g. 0.5, 4e-4, 1.2e3)
 FLOAT_RE = r"[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?"
 
-_TSIM_GATES = {"R_X", "R_Y", "R_Z", "U3"}
+_TSIM_GATES = {"R_X", "R_Y", "R_Z", "R_XX", "R_YY", "R_ZZ", "R_PAULI", "U3"}
 _GATE_NOT_FOUND_RE = re.compile(r"Gate not found: '(\w+)'")
-_GATE_USAGE_RE = re.compile(r"(?<!\[)\b(R_[A-Z]\([^)]*\)|R_[XYZ]\b|U3\([^)]*\)|U3\b)")
+_GATE_USAGE_RE = re.compile(
+    r"(?<!\[)\b(R_(?:XX|YY|ZZ|PAULI|[XYZ])\([^)]*\)|R_(?:XX|YY|ZZ|PAULI|[XYZ])\b|U3\([^)]*\)|U3\b)"
+)
 
 
 def enriched_stim_error(exc: ValueError, converted_text: str) -> ValueError:
@@ -38,6 +40,10 @@ def shorthand_to_stim(text: str) -> str:
         R_Z(0.3) 0          → I[R_Z(theta=0.3*pi)] 0
         R_X(0.25) 0         → I[R_X(theta=0.25*pi)] 0
         R_Y(-0.5) 0         → I[R_Y(theta=-0.5*pi)] 0
+        R_XX(0.25) 0 1      → SPP[R_XX(theta=0.25*pi)] X0*X1
+        R_YY(0.25) 0 1      → SPP[R_YY(theta=0.25*pi)] Y0*Y1
+        R_ZZ(0.25) 0 1      → SPP[R_ZZ(theta=0.25*pi)] Z0*Z1
+        R_PAULI(0.25) X0*Y1 → SPP[R_PAULI(theta=0.25*pi)] X0*Y1
         U3(0.3, 0.24, 0.49) 0 → I[U3(theta=0.3*pi, phi=0.24*pi, lambda=0.49*pi)] 0
 
     """
@@ -48,6 +54,36 @@ def shorthand_to_stim(text: str) -> str:
     text = re.sub(r"(?<!\[)\bT_DAG\b(?!\[)", "S_DAG[T]", text)
     text = re.sub(r"(?<!\[)\bT\b(?!\[)", "S[T]", text)
 
+    # R_XX/R_YY/R_ZZ must come before R_X/R_Y/R_Z to avoid partial matches
+    def replace_r_pp(m: re.Match) -> str:
+        pauli = m.group(1)  # "XX", "YY", or "ZZ"
+        angle = float(m.group(2))
+        q0 = m.group(3)
+        q1 = m.group(4)
+        if q0 == q1:
+            raise ValueError(
+                f"Duplicate target qubits in R_{pauli}: both targets are qubit {q0}."
+            )
+        p = pauli[0]
+        return f"SPP[R_{pauli}(theta={angle}*pi)] {p}{q0}*{p}{q1}"
+
+    text = re.sub(
+        rf"\bR_(XX|YY|ZZ)\(({FLOAT_RE})\)\s+(\d+)\s+(\d+)",
+        replace_r_pp,
+        text,
+    )
+
+    def replace_r_pauli(m: re.Match) -> str:
+        angle = float(m.group(1))
+        targets = m.group(2)
+        return f"SPP[R_PAULI(theta={angle}*pi)] {targets}"
+
+    text = re.sub(
+        rf"\bR_PAULI\(({FLOAT_RE})\)\s+((?:[XYZ]\d+)(?:\*[XYZ]\d+)*)",
+        replace_r_pauli,
+        text,
+    )
+
     def replace_rotation(m: re.Match) -> str:
         axis = m.group(1)
         return f"I[R_{axis}(theta={float(m.group(2))}*pi)]"
@@ -86,6 +122,10 @@ def stim_to_shorthand(text: str) -> str:
     Rewrites:
     - I[U3(theta=θ*pi, phi=φ*pi, lambda=λ*pi)] → U3(θ, φ, λ)
     - I[R_X(theta=α*pi)] / I[R_Y(...)] / I[R_Z(...)] → R_X(α) / R_Y(α) / R_Z(α)
+    - SPP[R_XX(theta=α*pi)] X0*X1 → R_XX(α) 0 1
+    - SPP[R_YY(theta=α*pi)] Y0*Y1 → R_YY(α) 0 1
+    - SPP[R_ZZ(theta=α*pi)] Z0*Z1 → R_ZZ(α) 0 1
+    - SPP[R_PAULI(theta=α*pi)] X0*Y1*Z2 → R_PAULI(α) X0*Y1*Z2
     - SPP[T] → TPP
     - SPP_DAG[T] → TPP_DAG
     - S[T] → T
@@ -115,6 +155,32 @@ def replace_rotation(m: re.Match) -> str:
         text,
     )
 
+    # Replace SPP[R_XX/R_YY/R_ZZ(...)] with R_XX/R_YY/R_ZZ(...)
+    def replace_spp_r_pp(m: re.Match) -> str:
+        pauli = m.group(1)
+        angle = m.group(2)
+        q0 = m.group(3)
+        q1 = m.group(4)
+        return f"R_{pauli}({angle}) {q0} {q1}"
+
+    text = re.sub(
+        rf"\bSPP\[R_(XX|YY|ZZ)\(theta=({FLOAT_RE})\*pi\)\]\s+[XYZ](\d+)\*[XYZ](\d+)",
+        replace_spp_r_pp,
+        text,
+    )
+
+    # Replace SPP[R_PAULI(...)] with R_PAULI(...)
+    def replace_spp_r_pauli(m: re.Match) -> str:
+        angle = m.group(1)
+        targets = m.group(2)
+        return f"R_PAULI({angle}) {targets}"
+
+    text = re.sub(
+        rf"\bSPP\[R_PAULI\(theta=({FLOAT_RE})\*pi\)\]\s+((?:[XYZ]\d+)(?:\*[XYZ]\d+)*)",
+        replace_spp_r_pauli,
+        text,
+    )
+
     # Replace SPP[T] and SPP_DAG[T] with TPP and TPP_DAG
     # Must come before S[T]/S_DAG[T] to avoid partial matches
     text = re.sub(r"(?<!\w)SPP_DAG\[T\](?!\w)", "TPP_DAG", text)
diff --git a/test/helpers/gate_matrices.py b/test/helpers/gate_matrices.py
index f669f196..15c48a9b 100644
--- a/test/helpers/gate_matrices.py
+++ b/test/helpers/gate_matrices.py
@@ -114,6 +114,31 @@
     "R_Z": lambda frac: np.array(
         [[np.exp(-1j * np.pi / 2 * frac), 0], [0, np.exp(1j * np.pi / 2 * frac)]]
     ),
+    # Two-qubit Pauli rotations: exp(-i * frac * pi/2 * PP)
+    "R_XX": lambda frac: np.array(
+        [
+            [np.cos(frac * np.pi / 2), 0, 0, -1j * np.sin(frac * np.pi / 2)],
+            [0, np.cos(frac * np.pi / 2), -1j * np.sin(frac * np.pi / 2), 0],
+            [0, -1j * np.sin(frac * np.pi / 2), np.cos(frac * np.pi / 2), 0],
+            [-1j * np.sin(frac * np.pi / 2), 0, 0, np.cos(frac * np.pi / 2)],
+        ]
+    ),
+    "R_YY": lambda frac: np.array(
+        [
+            [np.cos(frac * np.pi / 2), 0, 0, 1j * np.sin(frac * np.pi / 2)],
+            [0, np.cos(frac * np.pi / 2), -1j * np.sin(frac * np.pi / 2), 0],
+            [0, -1j * np.sin(frac * np.pi / 2), np.cos(frac * np.pi / 2), 0],
+            [1j * np.sin(frac * np.pi / 2), 0, 0, np.cos(frac * np.pi / 2)],
+        ]
+    ),
+    "R_ZZ": lambda frac: np.array(
+        [
+            [np.exp(-1j * frac * np.pi / 2), 0, 0, 0],
+            [0, np.exp(1j * frac * np.pi / 2), 0, 0],
+            [0, 0, np.exp(1j * frac * np.pi / 2), 0],
+            [0, 0, 0, np.exp(-1j * frac * np.pi / 2)],
+        ]
+    ),
     "U3": lambda frac_theta, frac_phi, frac_lambda: np.array(
         [
             [
diff --git a/test/integration/test_gate_unitaries.py b/test/integration/test_gate_unitaries.py
index 2e14fe54..cbda58d2 100644
--- a/test/integration/test_gate_unitaries.py
+++ b/test/integration/test_gate_unitaries.py
@@ -262,6 +262,43 @@ def test_tpp_z_equals_t():
     assert np.allclose(unitary, t_matrix)
 
 
+@pytest.mark.parametrize("instruction", ["R_XX", "R_YY", "R_ZZ"])
+def test_two_qubit_rot_instructions(instruction: str):
+    c = Circuit(f"""
+        R 0 1 2 3
+        H 0 1
+        CNOT 0 2 1 3
+        {instruction}(0.345) 0 1
+        M 0 1 2 3
+        """)
+    sampler = CompiledStateProbs(c)
+    mat = get_matrix(sampler)
+    expected = np.abs(ROT_GATE_MATRICES[instruction](0.345)) ** 2
+    assert np.allclose(mat, expected)
+
+
+@pytest.mark.parametrize(
+    "pauli_str, ref_key",
+    [
+        ("X0*X1", "R_XX"),
+        ("Y0*Y1", "R_YY"),
+        ("Z0*Z1", "R_ZZ"),
+    ],
+)
+def test_r_pauli_matches_two_qubit(pauli_str: str, ref_key: str):
+    c = Circuit(f"""
+        R 0 1 2 3
+        H 0 1
+        CNOT 0 2 1 3
+        R_PAULI(0.345) {pauli_str}
+        M 0 1 2 3
+        """)
+    sampler = CompiledStateProbs(c)
+    mat = get_matrix(sampler)
+    expected = np.abs(ROT_GATE_MATRICES[ref_key](0.345)) ** 2
+    assert np.allclose(mat, expected)
+
+
 def test_u3_instruction():
     c = Circuit("""
         R 0 1
diff --git a/test/unit/core/test_instructions.py b/test/unit/core/test_instructions.py
index 1bf214f5..011a7e18 100644
--- a/test/unit/core/test_instructions.py
+++ b/test/unit/core/test_instructions.py
@@ -4,7 +4,7 @@
     SINGLE_QUBIT_GATE_MATRICES,
     TWO_QUBIT_GATE_MATRICES,
 )
-from typing import Callable
+from typing import Callable, Literal
 
 import numpy as np
 import pytest
@@ -98,6 +98,21 @@ def test_u3(frac_theta: Fraction, frac_phi: Fraction, frac_lambda: Fraction):
     assert np.allclose(result, expected)
 
 
+@pytest.mark.parametrize("frac", [Fraction(1, 5), Fraction(-1, 3), Fraction(1, 7)])
+@pytest.mark.parametrize("pauli_char", ["X", "Y", "Z"])
+def test_r_pauli_two_qubit(frac: Fraction, pauli_char: Literal["X", "Y", "Z"]):
+    paulis: list[tuple[Literal["X", "Y", "Z"], int]] = [
+        (pauli_char, 0),
+        (pauli_char, 1),
+    ]
+    b = GraphRepresentation()
+    instructions.r_pauli(b, paulis, frac)
+    b.graph.normalize()
+    result = b.graph.to_matrix()
+    expected = ROT_GATE_MATRICES[f"R_{pauli_char}{pauli_char}"](frac)
+    assert np.allclose(result, expected)
+
+
 @pytest.mark.parametrize(
     "gate_func, matrix",
     [
diff --git a/test/unit/core/test_parse.py b/test/unit/core/test_parse.py
index 1d838603..198be241 100644
--- a/test/unit/core/test_parse.py
+++ b/test/unit/core/test_parse.py
@@ -975,6 +975,13 @@ def test_negative_angle_parsed(self):
             {"theta": Fraction(-1, 2)},
         )
 
+    def test_r_pp_tag_returns_theta(self):
+        for axis in ("R_XX", "R_YY", "R_ZZ", "R_PAULI"):
+            assert parse_parametric_tag(_instr(f"{axis}(theta=0.25*pi)")) == (
+                axis,
+                {"theta": Fraction(1, 4)},
+            )
+
     def test_non_parametric_tag_returns_none(self):
         # Tags without the name(...) shape are not parametric-looking; these are
         # used for non-parametric annotations like S[T] / SPP[T].
diff --git a/test/unit/test_circuit.py b/test/unit/test_circuit.py
index 061e3563..f7d42ca1 100644
--- a/test/unit/test_circuit.py
+++ b/test/unit/test_circuit.py
@@ -850,6 +850,95 @@ def test_append():
     assert "U3(0.3, 0.24, 0.49) 0" in str(c)
 
 
+def test_append_r_xx():
+    c = Circuit()
+    c.append("R_XX", [0, 1], arg=0.25)
+    assert "R_XX(0.25) 0 1" in str(c)
+
+
+def test_append_r_yy():
+    c = Circuit()
+    c.append("R_YY", [2, 3], arg=-0.5)
+    assert "R_YY(-0.5) 2 3" in str(c)
+
+
+def test_append_r_zz():
+    c = Circuit()
+    c.append("R_ZZ", [0, 1], arg=[0.3])
+    assert "R_ZZ(0.3) 0 1" in str(c)
+
+
+def test_append_r_pauli():
+    c = Circuit()
+    c.append(
+        "R_PAULI",
+        [stim.target_x(0), stim.target_combiner(), stim.target_y(1)],
+        arg=0.25,
+    )
+    assert "R_PAULI(0.25) X0*Y1" in str(c)
+
+
+def test_append_r_xx_duplicate_qubits_raises():
+    c = Circuit()
+    with pytest.raises(ValueError, match="Duplicate target qubits"):
+        c.append("R_XX", [3, 3], arg=0.5)
+
+
+def test_append_r_xx_wrong_target_count_raises():
+    c = Circuit()
+    with pytest.raises(ValueError, match="exactly two"):
+        c.append("R_XX", [0], arg=0.5)
+
+
+def test_r_xx_shorthand():
+    c1 = Circuit("R_XX(0.25) 0 1")
+    c2 = Circuit("SPP[R_XX(theta=0.25*pi)] X0*X1")
+    assert c1._stim_circ == c2._stim_circ
+
+
+def test_r_pauli_shorthand():
+    c1 = Circuit("R_PAULI(0.25) X0*Y1")
+    c2 = Circuit("SPP[R_PAULI(theta=0.25*pi)] X0*Y1")
+    assert c1._stim_circ == c2._stim_circ
+
+
+def test_inverse_r_xx():
+    c = Circuit("R_XX(0.25) 0 1")
+    c_inv = c.inverse()
+    combined = (c + c_inv).to_matrix()
+    assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
+
+
+def test_inverse_r_yy():
+    c = Circuit("R_YY(-0.3) 0 1")
+    c_inv = c.inverse()
+    combined = (c + c_inv).to_matrix()
+    assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
+
+
+def test_inverse_r_zz():
+    c = Circuit("R_ZZ(0.5) 0 1")
+    c_inv = c.inverse()
+    combined = (c + c_inv).to_matrix()
+    assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
+
+
+def test_inverse_r_pauli():
+    c = Circuit("R_PAULI(0.25) X0*Y1*Z2")
+    c_inv = c.inverse()
+    combined = (c + c_inv).to_matrix()
+    assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
+
+
+def test_inverse_spp_dag_parametric():
+    c = Circuit.from_stim_program(
+        stim.Circuit("SPP_DAG[R_PAULI(theta=0.25*pi)] X0")
+    )
+    c_inv = c.inverse()
+    combined = (c + c_inv).to_matrix()
+    assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
+
+
 def test_append_circuit_instruction():
     c = Circuit()
     c.append(stim.CircuitInstruction("H", [0]))
diff --git a/test/unit/utils/test_clifford.py b/test/unit/utils/test_clifford.py
index f5d4c1b7..fbff5438 100644
--- a/test/unit/utils/test_clifford.py
+++ b/test/unit/utils/test_clifford.py
@@ -185,6 +185,23 @@ def test_multi_target_rotation(self):
         assert "I[" not in stim_str
         assert "Z 0 1 2" in stim_str
 
+    def test_spp_dag_parametric_expanded(self):
+        c = Circuit.from_stim_program(
+            stim.Circuit("SPP_DAG[R_PAULI(theta=0.5*pi)] X0*Y1")
+        )
+        stim_str = str(c.stim_circuit)
+        assert "SPP_DAG X0*Y1" in stim_str
+        assert "SPP X0*Y1" not in stim_str
+
+    def test_spp_dag_parametric_matches_parsed_unitary(self):
+        c = Circuit.from_stim_program(
+            stim.Circuit("SPP_DAG[R_PAULI(theta=0.5*pi)] X0*Y1")
+        )
+        expected = Circuit.from_stim_program(stim.Circuit("SPP_DAG X0*Y1"))
+        assert _unitaries_equal_up_to_global_phase(
+            c.to_matrix(), expected.to_matrix()
+        )
+
 
 class TestIsCliffordUserTags:
     """Identity instructions with non-parametric user tags (e.g. ``I[mytag]``)
diff --git a/test/unit/utils/test_program_text.py b/test/unit/utils/test_program_text.py
index 14a8dfbb..8117ec09 100644
--- a/test/unit/utils/test_program_text.py
+++ b/test/unit/utils/test_program_text.py
@@ -68,6 +68,57 @@ def test_shorthand_roundtrip():
     assert stim_to_shorthand(shorthand_to_stim(text)) == text
 
 
+def test_shorthand_to_stim_r_xx():
+    text = "R_XX(0.25) 0 1"
+    expected = "SPP[R_XX(theta=0.25*pi)] X0*X1"
+    assert shorthand_to_stim(text) == expected
+
+
+def test_shorthand_to_stim_r_yy():
+    text = "R_YY(-0.5) 2 3"
+    expected = "SPP[R_YY(theta=-0.5*pi)] Y2*Y3"
+    assert shorthand_to_stim(text) == expected
+
+
+def test_shorthand_to_stim_r_zz():
+    text = "R_ZZ(0.3) 0 1"
+    expected = "SPP[R_ZZ(theta=0.3*pi)] Z0*Z1"
+    assert shorthand_to_stim(text) == expected
+
+
+def test_shorthand_to_stim_r_pauli():
+    text = "R_PAULI(0.25) X0*Y1*Z2"
+    expected = "SPP[R_PAULI(theta=0.25*pi)] X0*Y1*Z2"
+    assert shorthand_to_stim(text) == expected
+
+
+def test_stim_to_shorthand_r_xx():
+    text = "SPP[R_XX(theta=0.25*pi)] X0*X1"
+    expected = "R_XX(0.25) 0 1"
+    assert stim_to_shorthand(text) == expected
+
+
+def test_stim_to_shorthand_r_pauli():
+    text = "SPP[R_PAULI(theta=0.25*pi)] X0*Y1*Z2"
+    expected = "R_PAULI(0.25) X0*Y1*Z2"
+    assert stim_to_shorthand(text) == expected
+
+
+def test_shorthand_roundtrip_r_xx():
+    text = "R_XX(0.25) 0 1"
+    assert stim_to_shorthand(shorthand_to_stim(text)) == text
+
+
+def test_shorthand_roundtrip_r_pauli():
+    text = "R_PAULI(0.25) X0*Y1*Z2"
+    assert stim_to_shorthand(shorthand_to_stim(text)) == text
+
+
+def test_shorthand_r_xx_duplicate_qubits_raises():
+    with pytest.raises(ValueError, match="Duplicate target qubits"):
+        shorthand_to_stim("R_XX(0.5) 3 3")
+
+
 def test_shorthand_scientific_notation():
     result = shorthand_to_stim("R_Z(4e-4) 0")
     assert "I[R_Z(theta=0.0004*pi)]" in result