Fix inconsistent Pauli gate multiplication

cirq-core/cirq/ops/pauli_string.py

@@ -271,6 +271,8 @@ def __mul__(self, other):
             known = True
         elif isinstance(other, (PauliString, numbers.Number)):
             known = True
+        elif (as_pauli_string := _try_interpret_as_pauli_string(other)) is not None:
+            return self * as_pauli_string
         if known:
             return PauliString(
                 cast(PAULI_STRING_LIKE, other),
@@ -297,6 +299,8 @@ def __rmul__(self, other) -> PauliString:
 
         if isinstance(other, raw_types.Operation) and isinstance(other.gate, identity.IdentityGate):
             return self  # pragma: no cover
+        elif (as_pauli_string := _try_interpret_as_pauli_string(other)) is not None:
+            return as_pauli_string * self
 
         # Note: PauliString case handled by __mul__.
         return NotImplemented

cirq-core/cirq/ops/raw_types.py

@@ -483,10 +483,62 @@ def _commutes_(self, other: Any, *, atol: float = 1e-8) -> None | NotImplemented
 
     def _mul_with_qubits(self, qubits: tuple[cirq.Qid, ...], other):
         """cirq.GateOperation.__mul__ delegates to this method."""
+        if isinstance(other, Operation):
+            try:
+                # Try using pauli expansion if both operations have single-item expansions
+                pauli_expansion_self = protocols.pauli_expansion(self.on(*qubits))
+                pauli_expansion_other = protocols.pauli_expansion(other)
+
+                if (
+                    pauli_expansion_self is not None
+                    and len(pauli_expansion_self) == 1
+                    and pauli_expansion_other is not None
+                    and len(pauli_expansion_other) == 1
+                ):
+
+                    gate_self, coef_self = next(iter(pauli_expansion_self.items()))
+                    gate_other, coef_other = next(iter(pauli_expansion_other.items()))
+
+                    from cirq.ops.pauli_string import PauliString
+
+                    return (
+                        coef_self
+                        * PauliString({q: gate_self for q in qubits})
+                        * coef_other
+                        * PauliString({q: gate_other for q in qubits})
+                    )
+            except TypeError:
+                return NotImplemented
         return NotImplemented
 
     def _rmul_with_qubits(self, qubits: tuple[cirq.Qid, ...], other):
         """cirq.GateOperation.__rmul__ delegates to this method."""
+        if isinstance(other, Operation):
+            try:
+                # Try using pauli expansion if both operations have single-item expansions
+                pauli_expansion_self = protocols.pauli_expansion(self.on(*qubits))
+                pauli_expansion_other = protocols.pauli_expansion(other)
+
+                if (
+                    pauli_expansion_self is not None
+                    and len(pauli_expansion_self) == 1
+                    and pauli_expansion_other is not None
+                    and len(pauli_expansion_other) == 1
+                ):
+
+                    gate_self, coef_self = next(iter(pauli_expansion_self.items()))
+                    gate_other, coef_other = next(iter(pauli_expansion_other.items()))
+
+                    from cirq.ops.pauli_string import PauliString
+
+                    return (
+                        coef_other
+                        * PauliString({q: gate_other for q in qubits})
+                        * coef_self
+                        * PauliString({q: gate_self for q in qubits})
+                    )
+            except TypeError:
+                return NotImplemented
         return NotImplemented
 
     def _json_dict_(self) -> dict[str, Any]: