Repository: sympy/sympy

Issue:
Quantum Operators * their inverse become scalar 1 instead of IdentityOperator
### Situation ###
Operator classes like `Matrix` or `MatrixSymbol` represent scalar results as multiple of their respective Identity operator, i.e. operator * operator -> operator (except for 0):
```
from sympy.matrices import *
M = MatrixSymbol('M', 3, 3)
print (M * M**-1, M * M.inv(), type(pow(M, 0)) 
#->  Identity(3)  Identity(3)  sympy.matrices.expressions.special.Identity
```
However operators from the `sympy.physics.quantum Operator` class return scalars, though they are also linear operators:
```
from sympy.physics.quantum import *
A = Operator('A')
print(A * A**-1, A * A.inv(),  type(pow(A, 0)))
#->  1   1    <class 'sympy.core.numbers.One'>
```
### Assessment ###
Mathematically speaking this is not relevant. 
However it is unusual - even in mathematical notation one typically writes A * A**-1 = I  - and it is inconsistent within SymPy: 
Since SymPy typically has `Identity.is_commutative==False`,  `args_cnc` keeps the Identity in the nc_part, so code can rely on that the product of two operators is again an operator type (except for 0). With quantum.Operators all expression handling code that directly or indirectly relies on args_cnc behaves differently, which is inconvenient. 
As another example today's output of _TensorProduct(A * A**-1, B * B.inv())_ -> _1 x 1_ is much less intuitive than the expected  _I x I_  .  
Furthermore: If it was decided that the result should be represented as _I x I_, the term _1 x 1_ would be reserved for _really_ scalar results like  _(<1| x <0| ) x (|1> x |0>) = 1 x 1 = 1_ .

### Proposition for Fix ###
Ensuring `pow(A, 0)==IdentityOperator()` gives consistent behaviour with mathematical usances and other Sympy packages.
I have checked this can be implemented with minimal effort by adding a `__pow__()` method to the `quantum.Operator`  class that just adds handling of the `exponent==0` case and can provide a PR if requested.


Discussion:

Patch:
diff --git a/doc/src/explanation/active-deprecations.md b/doc/src/explanation/active-deprecations.md
index ec8ad588a3..d8b1c469fb 100644
--- a/doc/src/explanation/active-deprecations.md
+++ b/doc/src/explanation/active-deprecations.md
@@ -76,7 +76,25 @@ SymPy deprecation warnings.
 
 ## Version 1.14
 
-There are no deprecations yet for SymPy 1.14.
+(deprecated-operator-identity)=
+### Deprecated IdentityOperator from physics.quantum
+
+The ``IdentityOperator`` in the ``sympy.physics.quantum`` moddule has been
+deprecated. Originally, we thought that it would be helpful to have a
+multiplicative identity for quantum operators and states. However, at this
+time, it is unused in `sympy.physics.quantum` for anything other than tests
+of its own behavior. In addition, users were finding inconsistencies in 
+the behavior of ``IdentityOperator`` compared to what is expected by a
+multiplicative identity.
+
+Moving forward, we recommend that users use the scalar `S.One` as the 
+multiplicative identity for all operators and states in the quantum
+module. The code in ``sympy.physics.quantum`` currently does not ever
+return an ``IdentityOperator`` so the only place users will encounter
+its usage is in their own code.
+
+The existing implementation will remain, along with its tests for at least
+one year after the 1.14 release.
 
 ## Version 1.13
 
diff --git a/sympy/physics/quantum/boson.py b/sympy/physics/quantum/boson.py
index 3be2ebc45c..4dfd2286b1 100644
--- a/sympy/physics/quantum/boson.py
+++ b/sympy/physics/quantum/boson.py
@@ -7,7 +7,7 @@
 from sympy.functions.elementary.exponential import exp
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.physics.quantum import Operator
-from sympy.physics.quantum import HilbertSpace, FockSpace, Ket, Bra, IdentityOperator
+from sympy.physics.quantum import HilbertSpace, FockSpace, Ket, Bra
 from sympy.functions.special.tensor_functions import KroneckerDelta
 
 
@@ -94,9 +94,6 @@ def _eval_adjoint(self):
 
     def __mul__(self, other):
 
-        if other == IdentityOperator(2):
-            return self
-
         if isinstance(other, Mul):
             args1 = tuple(arg for arg in other.args if arg.is_commutative)
             args2 = tuple(arg for arg in other.args if not arg.is_commutative)
diff --git a/sympy/physics/quantum/operator.py b/sympy/physics/quantum/operator.py
index 8c540dc016..d5869a1607 100644
--- a/sympy/physics/quantum/operator.py
+++ b/sympy/physics/quantum/operator.py
@@ -20,6 +20,7 @@
 from sympy.physics.quantum.dagger import Dagger
 from sympy.physics.quantum.qexpr import QExpr, dispatch_method
 from sympy.matrices import eye
+from sympy.utilities.exceptions import sympy_deprecation_warning
 
 __all__ = [
     'Operator',
@@ -258,6 +259,10 @@ class IdentityOperator(Operator):
     """An identity operator I that satisfies op * I == I * op == op for any
     operator op.
 
+    .. deprecated:: 1.14.
+        Use the scalar S.One instead as the multiplicative identity for
+        operators and states.
+
     Parameters
     ==========
 
@@ -283,6 +288,14 @@ def default_args(self):
         return (oo,)
 
     def __init__(self, *args, **hints):
+        sympy_deprecation_warning(
+            """
+            IdentityOperator has been deprecated. In the future, please use
+            S.One as the identity for quantum operators and states.
+            """,
+            deprecated_since_version="1.14",
+            active_deprecations_target='deprecated-operator-identity',
+        )
         if not len(args) in (0, 1):
             raise ValueError('0 or 1 parameters expected, got %s' % args)
 


Test patch:
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index 8969b6d1ec..b617f414de 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -18,7 +18,7 @@
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import sin
 
-from sympy.testing.pytest import SKIP
+from sympy.testing.pytest import SKIP, warns_deprecated_sympy
 
 a, b, c, x, y, z, s = symbols('a,b,c,x,y,z,s')
 
@@ -3754,8 +3754,9 @@ def test_sympy__physics__quantum__operator__HermitianOperator():
 
 
 def test_sympy__physics__quantum__operator__IdentityOperator():
-    from sympy.physics.quantum.operator import IdentityOperator
-    assert _test_args(IdentityOperator(5))
+    with warns_deprecated_sympy():
+        from sympy.physics.quantum.operator import IdentityOperator
+        assert _test_args(IdentityOperator(5))
 
 
 def test_sympy__physics__quantum__operator__Operator():
diff --git a/sympy/physics/quantum/tests/test_dagger.py b/sympy/physics/quantum/tests/test_dagger.py
index 0d379095de..8aff96506d 100644
--- a/sympy/physics/quantum/tests/test_dagger.py
+++ b/sympy/physics/quantum/tests/test_dagger.py
@@ -7,7 +7,7 @@
 
 from sympy.physics.quantum.dagger import adjoint, Dagger
 from sympy.external import import_module
-from sympy.testing.pytest import skip
+from sympy.testing.pytest import skip, warns_deprecated_sympy
 from sympy.physics.quantum.operator import Operator, IdentityOperator
 
 
@@ -37,9 +37,10 @@ def test_matrix():
 
 def test_dagger_mul():
     O = Operator('O')
-    I = IdentityOperator()
     assert Dagger(O)*O == Dagger(O)*O
-    assert Dagger(O)*O*I == Mul(Dagger(O), O)*I
+    with warns_deprecated_sympy():
+        I = IdentityOperator()
+        assert Dagger(O)*O*I == Mul(Dagger(O), O)*I
     assert Dagger(O)*Dagger(O) == Dagger(O)**2
     assert Dagger(O)*Dagger(I) == Dagger(O)
 
diff --git a/sympy/physics/quantum/tests/test_operator.py b/sympy/physics/quantum/tests/test_operator.py
index df281fa833..8950fc9b93 100644
--- a/sympy/physics/quantum/tests/test_operator.py
+++ b/sympy/physics/quantum/tests/test_operator.py
@@ -17,6 +17,8 @@
 from sympy.physics.quantum.trace import Tr
 from sympy.matrices import eye
 
+from sympy.testing.pytest import warns_deprecated_sympy
+
 
 class CustomKet(Ket):
     @classmethod
@@ -53,7 +55,8 @@ def test_operator():
     assert t_op.label[0] == Symbol(t_op.default_args()[0])
 
     assert Operator() == Operator("O")
-    assert A*IdentityOperator() == A
+    with warns_deprecated_sympy():
+        assert A*IdentityOperator() == A
 
 
 def test_operator_inv():
@@ -88,27 +91,28 @@ def test_unitary():
 
 
 def test_identity():
-    I = IdentityOperator()
-    O = Operator('O')
-    x = Symbol("x")
-
-    assert isinstance(I, IdentityOperator)
-    assert isinstance(I, Operator)
-
-    assert I * O == O
-    assert O * I == O
-    assert I * Dagger(O) == Dagger(O)
-    assert Dagger(O) * I == Dagger(O)
-    assert isinstance(I * I, IdentityOperator)
-    assert isinstance(3 * I, Mul)
-    assert isinstance(I * x, Mul)
-    assert I.inv() == I
-    assert Dagger(I) == I
-    assert qapply(I * O) == O
-    assert qapply(O * I) == O
-
-    for n in [2, 3, 5]:
-        assert represent(IdentityOperator(n)) == eye(n)
+    with warns_deprecated_sympy():
+        I = IdentityOperator()
+        O = Operator('O')
+        x = Symbol("x")
+
+        assert isinstance(I, IdentityOperator)
+        assert isinstance(I, Operator)
+
+        assert I * O == O
+        assert O * I == O
+        assert I * Dagger(O) == Dagger(O)
+        assert Dagger(O) * I == Dagger(O)
+        assert isinstance(I * I, IdentityOperator)
+        assert isinstance(3 * I, Mul)
+        assert isinstance(I * x, Mul)
+        assert I.inv() == I
+        assert Dagger(I) == I
+        assert qapply(I * O) == O
+        assert qapply(O * I) == O
+
+        for n in [2, 3, 5]:
+            assert represent(IdentityOperator(n)) == eye(n)
 
 
 def test_outer_product():
diff --git a/sympy/physics/quantum/tests/test_qapply.py b/sympy/physics/quantum/tests/test_qapply.py
index b2d371570b..8394778224 100644
--- a/sympy/physics/quantum/tests/test_qapply.py
+++ b/sympy/physics/quantum/tests/test_qapply.py
@@ -17,6 +17,7 @@
 from sympy.physics.quantum.density import Density
 from sympy.physics.quantum.qubit import Qubit, QubitBra
 from sympy.physics.quantum.boson import BosonOp, BosonFockKet, BosonFockBra
+from sympy.testing.pytest import warns_deprecated_sympy
 
 
 j, jp, m, mp = symbols("j j' m m'")
@@ -146,5 +147,6 @@ def test_issue24158_ket_times_op():
     P2 = qapply(P2, dagger = True) # unsatisfactorily -> <0|*X(0), expect <1| since dagger=True
     assert qapply(P2, dagger = True) == QubitBra(1) # qapply(P1) -> 0 before fix
     # Pull Request 24237: IdentityOperator from the right without dagger=True option
-    assert qapply(QubitBra(1)*IdentityOperator()) == QubitBra(1)
-    assert qapply(IdentityGate(0)*(Qubit(0) + Qubit(1))) == Qubit(0) + Qubit(1)
+    with warns_deprecated_sympy():
+        assert qapply(QubitBra(1)*IdentityOperator()) == QubitBra(1)
+        assert qapply(IdentityGate(0)*(Qubit(0) + Qubit(1))) == Qubit(0) + Qubit(1)