Summation ranges for QM: i<=j versus i<j
**Description**
[class QuadraticModel](https://github.com/dwavesystems/dimod/blob/655e330d5813546630191a7b1e594e63d240a33e/dimod/quadratic/quadratic_model.py#L127) defines the summation range for quadratic coefficients as :math:`\\sum_{i<j} b_{i, j} x_i x_j`, and that is accurate for binary variables but quadratic models allow self loops for int variables, so it seems to me the summation for QMs should be :math:`\\sum_{i \\le j} b_{i, j} x_i x_j`

**Steps To Reproduce**

```
>>> from dimod import QuadraticModel
>>> qm = QuadraticModel({}, {(0, 0): 1}, 0, {0: "INTEGER"})
```

**Expected Behavior**
Update definition as needed.

**Additional Context**
https://docs.ocean.dwavesys.com/en/stable/concepts/cqm.html
diff --git a/dimod/constrained.py b/dimod/constrained.py
index 812b367a..2bfd51a1 100644
--- a/dimod/constrained.py
+++ b/dimod/constrained.py
@@ -72,19 +72,24 @@ class ConstrainedQuadraticModel:
 
         \\begin{align}
             \\text{Minimize an objective:} & \\\\
-            & \\sum_{i} a_i x_i + \\sum_{i<j} b_{ij} x_i x_j + c, \\\\
+            & \\sum_{i} a_i x_i + \\sum_{i \\le j} b_{ij} x_i x_j + c, \\\\
             \\text{Subject to constraints:} & \\\\
-            & \\sum_i a_i^{(c)} x_i + \\sum_{i<j} b_{ij}^{(c)} x_i x_j+ c^{(c)} \\le 0,
+            & \\sum_i a_i^{(c)} x_i + \\sum_{i \\le j} b_{ij}^{(c)} x_i x_j+ c^{(c)} \\le 0,
             \\quad c=1, \\dots, C_{\\rm ineq.}, \\\\
-            & \\sum_i a_i^{(d)} x_i + \\sum_{i<j} b_{ij}^{(d)} x_i x_j + c^{(d)} = 0,
+            & \\sum_i a_i^{(d)} x_i + \\sum_{i \\le j} b_{ij}^{(d)} x_i x_j + c^{(d)} = 0,
             \\quad d=1, \\dots, C_{\\rm eq.},
         \\end{align}
 
-    where :math:`\\{ x_i\\}_{i=1, \\dots, N}` can be binary or integer
+    where :math:`\\{ x_i\\}_{i=1, \\dots, N}` can be binary\\ [#]_ or integer
     variables, :math:`a_{i}, b_{ij}, c` are real values and
     :math:`C_{\\rm ineq.}, C_{\\rm eq,}` are the number of inequality and
     equality constraints respectively.
 
+    .. [#]
+        For binary variables, the range of the quadratic-term summation is
+        :math:`i < j` because :math:`x^2 = x` for binary values :math:`\\{0, 1\\}`
+        and :math:`s^2 = 1` for spin values :math:`\\{-1, 1\\}`.
+
     The objective and constraints are encoded as either :class:`.QuadraticModel`
     or :class:`.BinaryQuadraticModel` depending on the variable types used.
 
diff --git a/dimod/quadratic/quadratic_model.py b/dimod/quadratic/quadratic_model.py
index bf83edb5..3d177fd2 100644
--- a/dimod/quadratic/quadratic_model.py
+++ b/dimod/quadratic/quadratic_model.py
@@ -124,11 +124,16 @@ class QuadraticModel(QuadraticViewsMixin):
 
     .. math::
 
-        E(x) = \\sum_i a_i x_i + \\sum_{i<j} b_{i, j} x_i x_j + c
+        E(x) = \\sum_i a_i x_i + \\sum_{i \\le j} b_{i, j} x_i x_j + c
 
-    where :math:`\\{ x_i\\}_{i=1, \\dots, N}` can be binary or integer
+    where :math:`\\{ x_i\\}_{i=1, \\dots, N}` can be binary\\ [#]_ or integer
     variables and :math:`a_{i}, b_{ij}, c` are real values.
 
+    .. [#]
+        For binary variables, the range of the quadratic-term summation is
+        :math:`i < j` because :math:`x^2 = x` for binary values :math:`\\{0, 1\\}`
+        and :math:`s^2 = 1` for spin values :math:`\\{-1, 1\\}`.
+
     """
     _DATA_CLASSES = {
         np.dtype(np.float32): cyQM_float32,
['dimod/quadratic/quadratic_model.py', 'dimod/constrained.py']