Repository: sympy/sympy

Issue:
【SymPy.physics.quantum】 Thermal state in quantum optics undefined
Hi,

Given that coherent state was defined within this package by BosonCoherentKet(), I wonder why thermal state for photons is not defined.

**Coherent state**:
![image](https://github.com/sympy/sympy/assets/158588356/2de6fa0e-4ecc-4a9b-ac71-583d418ca2bc)
**Thermal state**:
![image](https://github.com/sympy/sympy/assets/158588356/516130d9-0550-4094-b1c9-02e5aedfafd8)
Although thermal state is a mixed state and should be defined in terms of density matrix, but I believe in principle this can be done.
I tried to define it by myself using BosonFockKet() and OuterProduct(), but unfortunately this doesn't work, because the Sum() function cannot recognize the photon number 'n' inside a ket as dummy index for summation.
# This method doesn't work, because the sum function doesn't recognize the dummy index inside the bracket notation.
`density_matrix_subspace = OuterProduct(ket_n,Dagger(ket_n))`
`density_matrix = Sum(nbar**n/(1+nbar)**(n+1)*density_matrix_subspace,(n, 0, oo)).doit()`

Any suggestion for this issue?
P.S. I think this is the same issue encountered in #23820: sympy can't handle symbolic Sum of JzKets.

Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index 5a0ac713f4..882fdc493e 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1174,7 +1174,6 @@ Prasoon Shukla <prasoon92.iitr@gmail.com>
 Prateek Papriwal <papriwalprateek@gmail.com>
 Praveen Perumal <ppraveen98841@gmail.com> Praveen Perumal <58477724+solvedbiscuit71@users.noreply.github.com>
 Praveen Sahu <povinsahu@gmail.com> povinsahu1909 <povinsahu@gmail.com>
-Prayag <v.prayag2005@gmail.com> vprayag2005 <v.prayag2005@gamil.com>
 Prayush Dawda <35144226+iamprayush@users.noreply.github.com>
 Prempal Singh <prempal.42@gmail.com>
 Prey Patel <patel.prey@iitgn.ac.in>
diff --git a/sympy/functions/elementary/integers.py b/sympy/functions/elementary/integers.py
index d40af92f7d..645bfc4bd7 100644
--- a/sympy/functions/elementary/integers.py
+++ b/sympy/functions/elementary/integers.py
@@ -219,7 +219,7 @@ def __ge__(self, other):
                 return self.args[0] >= other
             if other.is_number and other.is_real:
                 return self.args[0] >= ceiling(other)
-        if self.args[0] == other and other.is_real and other.is_noninteger:
+        if self.args[0] == other and other.is_real:
             return S.false
         if other is S.NegativeInfinity and self.is_finite:
             return S.true
@@ -247,8 +247,8 @@ def __lt__(self, other):
                 return self.args[0] < other
             if other.is_number and other.is_real:
                 return self.args[0] < ceiling(other)
-        if self.args[0] == other and other.is_real and other.is_noninteger:
-            return S.true
+        if self.args[0] == other and other.is_real:
+            return S.false
         if other is S.Infinity and self.is_finite:
             return S.true
 
@@ -387,8 +387,8 @@ def __gt__(self, other):
                 return self.args[0] > other
             if other.is_number and other.is_real:
                 return self.args[0] > floor(other)
-        if self.args[0] == other and other.is_real and other.is_noninteger:
-            return S.true
+        if self.args[0] == other and other.is_real:
+            return S.false
         if other is S.NegativeInfinity and self.is_finite:
             return S.true
 
@@ -415,7 +415,7 @@ def __le__(self, other):
                 return self.args[0] <= other
             if other.is_number and other.is_real:
                 return self.args[0] <= floor(other)
-        if self.args[0] == other and other.is_real and other.is_noninteger:
+        if self.args[0] == other and other.is_real:
             return S.false
         if other is S.Infinity and self.is_finite:
             return S.true
diff --git a/sympy/geometry/ellipse.py b/sympy/geometry/ellipse.py
index 199db25fde..3b100d91c6 100644
--- a/sympy/geometry/ellipse.py
+++ b/sympy/geometry/ellipse.py
@@ -1475,7 +1475,7 @@ def section_modulus(self, point=None):
 
 
 class Circle(Ellipse):
-    r"""A circle in space.
+    """A circle in space.
 
     Constructed simply from a center and a radius, from three
     non-collinear points, or the equation of a circle.
@@ -1526,7 +1526,7 @@ class Circle(Ellipse):
     (sqrt(2)/2, sqrt(2)/2, sqrt(2)/2, Point2D(1/2, 1/2))
 
     A circle can be constructed from an equation in the form
-    `ax^2 + by^2 + gx + hy + c = 0`, too:
+    `a*x**2 + by**2 + gx + hy + c = 0`, too:
 
     >>> Circle(x**2 + y**2 - 25)
     Circle(Point2D(0, 0), 5)
diff --git a/sympy/physics/quantum/qapply.py b/sympy/physics/quantum/qapply.py
index 2109ed1abc..87379c7e3e 100644
--- a/sympy/physics/quantum/qapply.py
+++ b/sympy/physics/quantum/qapply.py
@@ -4,6 +4,7 @@
 * Sometimes the final result needs to be expanded, we should do this by hand.
 """
 
+from sympy.concrete import Sum
 from sympy.core.add import Add
 from sympy.core.mul import Mul
 from sympy.core.power import Pow
@@ -46,6 +47,9 @@ def qapply(e, **options):
           (default: False).
         * ``ip_doit``: call ``.doit()`` in inner products when they are
           encountered (default: True).
+        * ``sum_doit``: call ``.doit()`` on sums when they are encountered
+          (default: False). This is helpful for collapsing sums over Kronecker
+          delta's that are created when calling ``qapply``.
 
     Returns
     =======
@@ -72,6 +76,7 @@ def qapply(e, **options):
     from sympy.physics.quantum.density import Density
 
     dagger = options.get('dagger', False)
+    sum_doit = options.get('sum_doit', False)
 
     if e == 0:
         return S.Zero
@@ -106,6 +111,13 @@ def qapply(e, **options):
     elif isinstance(e, TensorProduct):
         return TensorProduct(*[qapply(t, **options) for t in e.args])
 
+    # For a Sum, call qapply on its function.
+    elif isinstance(e, Sum):
+        result = Sum(qapply(e.function, **options), *e.limits)
+        if sum_doit:
+            result = result.doit()
+        return result
+
     # For a Pow, call qapply on its base.
     elif isinstance(e, Pow):
         return qapply(e.base, **options)**e.exp
@@ -133,6 +145,7 @@ def qapply(e, **options):
 def qapply_Mul(e, **options):
 
     ip_doit = options.get('ip_doit', True)
+    sum_doit = options.get('sum_doit', False)
 
     args = list(e.args)
 
@@ -177,10 +190,36 @@ def qapply_Mul(e, **options):
         result = TensorProduct(*[qapply(lhs.args[n]*rhs.args[n], **options) for n in range(len(lhs.args))]).expand(tensorproduct=True)
         return qapply_Mul(e.func(*args), **options)*result
 
+    # For Sums, move the Sum to the right.
+    if isinstance(rhs, Sum):
+        if isinstance(lhs, Sum):
+            if set(lhs.variables).intersection(set(rhs.variables)):
+                raise ValueError('Duplicated dummy indices in separate sums in qapply.')
+            limits = lhs.limits + rhs.limits
+            result = Sum(qapply(lhs.function*rhs.function, **options), *limits)
+            if sum_doit:
+                result = result.doit()
+            return qapply_Mul(e.func(*args)*result, **options)
+        else:
+            result = Sum(qapply(lhs*rhs.function, **options), rhs.limits)
+            if sum_doit:
+                result = result.doit()
+            return qapply_Mul(e.func(*args)*result, **options)
+
+    if isinstance(lhs, Sum):
+        result = Sum(qapply(lhs.function*rhs, **options), lhs.limits)
+        if sum_doit:
+            result = result.doit()
+        return qapply_Mul(e.func(*args)*result, **options)
+
     # Now try to actually apply the operator and build an inner product.
-    try:
-        result = lhs._apply_operator(rhs, **options)
-    except NotImplementedError:
+    _apply = getattr(lhs, '_apply_operator', None)
+    if _apply is not None:
+        try:
+            result = _apply(rhs, **options)
+        except NotImplementedError:
+            result = None
+    else:
         result = None
 
     if result is None:
diff --git a/sympy/physics/quantum/qexpr.py b/sympy/physics/quantum/qexpr.py
index 13f7f70294..64f7e2a200 100644
--- a/sympy/physics/quantum/qexpr.py
+++ b/sympy/physics/quantum/qexpr.py
@@ -97,10 +97,6 @@ class QExpr(Expr):
     # The separator used in printing the label.
     _label_separator = ''
 
-    @property
-    def free_symbols(self):
-        return {self}
-
     def __new__(cls, *args, **kwargs):
         """Construct a new quantum object.
 
diff --git a/sympy/physics/quantum/state.py b/sympy/physics/quantum/state.py
index 3688a54b4f..b2babef2e9 100644
--- a/sympy/physics/quantum/state.py
+++ b/sympy/physics/quantum/state.py
@@ -826,10 +826,6 @@ def _eval_conjugate(self):
     def _eval_transpose(self):
         return self
 
-    @property
-    def free_symbols(self):
-        return self.expr.free_symbols
-
     @property
     def is_commutative(self):
         """


Test patch:
diff --git a/sympy/functions/elementary/tests/test_integers.py b/sympy/functions/elementary/tests/test_integers.py
index a250e2f77a..1ac681a678 100644
--- a/sympy/functions/elementary/tests/test_integers.py
+++ b/sympy/functions/elementary/tests/test_integers.py
@@ -21,7 +21,6 @@
 i = Symbol('i', imaginary=True)
 y = Symbol('y', real=True)
 k, n = symbols('k,n', integer=True)
-b = Symbol('b', real=True, noninteger=True)
 
 
 def test_floor():
@@ -122,10 +121,10 @@ def test_floor():
     assert floor(factorial(50)/exp(1)) == \
         11188719610782480504630258070757734324011354208865721592720336800
 
-    assert (floor(y) < y).is_Relational
+    assert (floor(y) < y) == False
     assert (floor(y) <= y) == True
     assert (floor(y) > y) == False
-    assert (floor(y) >= y).is_Relational
+    assert (floor(y) >= y) == False
     assert (floor(x) <= x).is_Relational  # x could be non-real
     assert (floor(x) > x).is_Relational
     assert (floor(x) <= y).is_Relational  # arg is not same as rhs
@@ -134,10 +133,6 @@ def test_floor():
     assert (floor(y) < oo) == True
     assert (floor(y) >= -oo) == True
     assert (floor(y) > -oo) == True
-    assert (floor(b) < b) == True
-    assert (floor(b) <= b) == True
-    assert (floor(b) > b) == False
-    assert (floor(b) >= b) == False
 
     assert floor(y).rewrite(frac) == y - frac(y)
     assert floor(y).rewrite(ceiling) == -ceiling(-y)
@@ -316,9 +311,9 @@ def test_ceiling():
         11188719610782480504630258070757734324011354208865721592720336801
 
     assert (ceiling(y) >= y) == True
-    assert (ceiling(y) > y).is_Relational
+    assert (ceiling(y) > y) == False
     assert (ceiling(y) < y) == False
-    assert (ceiling(y) <= y).is_Relational
+    assert (ceiling(y) <= y) == False
     assert (ceiling(x) >= x).is_Relational  # x could be non-real
     assert (ceiling(x) < x).is_Relational
     assert (ceiling(x) >= y).is_Relational  # arg is not same as rhs
@@ -327,10 +322,6 @@ def test_ceiling():
     assert (ceiling(y) > -oo) == True
     assert (ceiling(y) <= oo) == True
     assert (ceiling(y) < oo) == True
-    assert (ceiling(b) < b) == False
-    assert (ceiling(b) <= b) == False
-    assert (ceiling(b) > b) == True
-    assert (ceiling(b) >= b) == True
 
     assert ceiling(y).rewrite(floor) == -floor(-y)
     assert ceiling(y).rewrite(frac) == y + frac(-y)
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 51b81775bb..7c7e950e04 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -2176,12 +2176,3 @@ def test_definite_integral_with_floats_issue_27231():
 
     # Verify that the result is approximately equal within a larger tolerance
     assert abs(result_numeric - expected_result.evalf()) < 1e-8
-
-
-def test_issue_27374():
-    #https://github.com/sympy/sympy/issues/27374
-    r = sqrt(x**2 + z**2)
-    u = erf(a*r/sqrt(2))/r
-    Ec = diff(u, z, z).subs([(x, sqrt(b*b-z*z))])
-    expected_result = -2*sqrt(2)*b*a**3*exp(-b**2*a**2/2)/(3*sqrt(pi))
-    assert simplify(integrate(Ec, (z, -b, b))) == expected_result
diff --git a/sympy/physics/quantum/tests/test_qexpr.py b/sympy/physics/quantum/tests/test_qexpr.py
index 667b826f68..c01817935a 100644
--- a/sympy/physics/quantum/tests/test_qexpr.py
+++ b/sympy/physics/quantum/tests/test_qexpr.py
@@ -1,11 +1,14 @@
 from sympy.core.numbers import Integer
 from sympy.core.symbol import Symbol
+from sympy.concrete import Sum
 from sympy.physics.quantum.qexpr import QExpr, _qsympify_sequence
 from sympy.physics.quantum.hilbert import HilbertSpace
 from sympy.core.containers import Tuple
 
 x = Symbol('x')
 y = Symbol('y')
+n = Symbol('n', integer=True)
+m = Symbol('m', integer=True)
 
 
 def test_qexpr_new():
@@ -32,12 +35,21 @@ def test_qexpr_commutative():
     q = QExpr._new_rawargs(Integer(0), Integer(1), HilbertSpace())
     assert q.is_commutative is False
 
-def test_qexpr_commutative_free_symbols():
-    q1 = QExpr(x)
-    assert q1.free_symbols.pop().is_commutative is False
 
-    q2 = QExpr('q2')
-    assert q2.free_symbols.pop().is_commutative is False
+def test_qexpr_free_symbols():
+    q1 = QExpr(x, y)
+    assert q1.free_symbols == {x, y}
+
+
+def test_qexpr_sum():
+    q1 = Sum(QExpr(n), (n,0,2))
+    assert q1.doit() == QExpr(0) + QExpr(1) + QExpr(2)
+
+    q2 = Sum(QExpr(n, m), (n, 0, 2), (m, 0, 2))
+    assert q2.doit() == QExpr(0, 0) + QExpr(0, 1) + QExpr(0, 2) + \
+        QExpr(1, 0) + QExpr(1, 1) + QExpr(1, 2) + \
+        QExpr(2, 0) + QExpr(2, 1) + QExpr(2, 2)
+
 
 def test_qexpr_subs():
     q1 = QExpr(x, y)
diff --git a/sympy/physics/quantum/tests/test_sho1d.py b/sympy/physics/quantum/tests/test_sho1d.py
index b9d9d37991..18d3862033 100644
--- a/sympy/physics/quantum/tests/test_sho1d.py
+++ b/sympy/physics/quantum/tests/test_sho1d.py
@@ -1,20 +1,27 @@
 """Tests for sho1d.py"""
 
+from sympy.concrete import Sum
+from sympy.core import oo
 from sympy.core.numbers import (I, Integer)
 from sympy.core.singleton import S
-from sympy.core.symbol import Symbol
+from sympy.core.symbol import Symbol, symbols
+from sympy.functions.combinatorial.factorials import factorial
+from sympy.functions.elementary.exponential import exp
 from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.functions.elementary.complexes import Abs
+from sympy.functions.special.tensor_functions import KroneckerDelta
 from sympy.physics.quantum import Dagger
 from sympy.physics.quantum.constants import hbar
 from sympy.physics.quantum import Commutator
 from sympy.physics.quantum.qapply import qapply
 from sympy.physics.quantum.innerproduct import InnerProduct
 from sympy.physics.quantum.cartesian import X, Px
-from sympy.functions.special.tensor_functions import KroneckerDelta
 from sympy.physics.quantum.hilbert import ComplexSpace
 from sympy.physics.quantum.represent import represent
+from sympy.simplify import simplify
 from sympy.external import import_module
-from sympy.testing.pytest import skip
+from sympy.tensor import IndexedBase, Idx
+from sympy.testing.pytest import skip, raises
 
 from sympy.physics.quantum.sho1d import (RaisingOp, LoweringOp,
                                         SHOKet, SHOBra,
@@ -33,6 +40,9 @@
 omega = Symbol('omega')
 m = Symbol('m')
 ndim = Integer(4)
+p = Symbol('p', is_integer=True)
+q = Symbol('q', nonnegative=True, is_integer=True)
+
 
 np = import_module('numpy')
 scipy = import_module('scipy', import_kwargs={'fromlist': ['sparse']})
@@ -118,3 +128,49 @@ def test_SHOKet():
     assert k.hilbert_space == ComplexSpace(S.Infinity)
     assert k3_rep[k3.n, 0] == Integer(1)
     assert b3_rep[0, b3.n] == Integer(1)
+
+def test_sho_sums():
+    e1 = Sum(SHOKet(p)*SHOBra(p), (p, 0, 1))
+    assert e1.doit() == SHOKet(0)*SHOBra(0) + SHOKet(1)*SHOBra(1)
+
+    # Test qapply with Sum on the left
+    assert qapply(
+        Sum(SHOKet(p)*SHOBra(p), (p, 0, oo))*SHOKet(q),
+        sum_doit=True
+    ) == SHOKet(q)
+
+    # Test qapply with Sum on the right
+    a = IndexedBase('a')
+    n = symbols('n', cls=Idx)
+    result = qapply(SHOBra(q)*Sum(a[n]*SHOKet(n), (n,0,oo)), sum_doit=True)
+    assert result == a[q]
+
+    # Test qapply with a product of Sums
+    result = qapply(
+        SHOBra(q)*Sum(SHOKet(p)*SHOBra(p), (p, 0, oo))*Sum(a[n]*SHOKet(n), (n,0,oo)),
+        sum_doit=True
+    )
+    assert result == a[q]
+
+    with raises(ValueError):
+        result = qapply(
+            SHOBra(q)*Sum(SHOKet(p)*SHOBra(p), (p, 0, oo))*Sum(a[p]*SHOKet(p), (p,0,oo)),
+            sum_doit=True
+        )
+
+def test_sho_coherant_state():
+    alpha = Symbol('alpha', is_complex=True)
+    cstate = exp(-Abs(alpha)**2/S(2))*Sum(((alpha**p)/sqrt(factorial(p)))*SHOKet(p), (p,0,oo))
+    # Projection onto the number eigenstate
+    assert qapply(SHOBra(q)*cstate, sum_doit=True) == exp(-Abs(alpha)**2/S(2))*alpha**q/sqrt(factorial(q))
+    # Ensure that the coherent state is an eigenstate of anihilation operator
+    assert simplify(qapply(SHOBra(q)*a*cstate, sum_doit=True)) == simplify(qapply(SHOBra(q)*alpha*cstate, sum_doit=True))
+
+def test_issue_26495():
+    nbar = Symbol('nbar', is_real=True, nonnegative=True)
+    n = Symbol('n', is_integer=True)
+    i = Symbol('i', is_integer=True, nonnegative=True)
+    j = Symbol('j', is_integer=True, nonnegative=True)
+    rho = (1/(1+nbar))*Sum((nbar/(1+nbar))**n*SHOKet(n)*SHOBra(n), (n,0,oo))
+    result = qapply(SHOBra(i)*rho*SHOKet(j), sum_doit=True)
+    assert simplify(result) == nbar**j*(nbar+1)**(-j-1)*KroneckerDelta(i,j)