diff --git a/gbasis/evals/electrostatic_potential.py b/gbasis/evals/electrostatic_potential.py
index d90dd3ca..5f3a90bd 100644
--- a/gbasis/evals/electrostatic_potential.py
+++ b/gbasis/evals/electrostatic_potential.py
@@ -14,6 +14,21 @@ def electrostatic_potential(
 ):
     r"""Return the electrostatic potentials of the basis set in the Cartesian form.
 
+    .. math::
+        \begin{align}
+        &- \left(
+            - \sum_A \frac{Z_A}{|\mathbf{R}_C - \mathbf{R}_A|}
+            + \sum_{ab} \gamma_{ab} \int \phi_a(\mathbf{r}) \frac{1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
+        \right)\\
+        &=
+        \sum_A \frac{Z_A}{|\mathbf{R}_C - \mathbf{R}_A|}
+        - \sum_{ab} \gamma_{ab} \int \phi_a(\mathbf{r}) \frac{1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
+        \end{align}
+
+    where :math:`\mathbf{R}_C` is the coordinate of a unitary point charge, :math:`\mathbf{R}_A` is the 
+    coordinate of the nucleus :math:`A`, :math:`Z_A` its charge, and :math:`\gamma_{ab}` is the 
+    one-electron density matrix.
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/angular_momentum.py b/gbasis/integrals/angular_momentum.py
index eb0c077c..f19046be 100644
--- a/gbasis/integrals/angular_momentum.py
+++ b/gbasis/integrals/angular_momentum.py
@@ -159,6 +159,21 @@ def construct_array_contraction(contractions_one, contractions_two):
 def angular_momentum_integral(basis, transform=None):
     r"""Return the integral over :math:`hat{L}` of the given basis set.
 
+    .. math::
+        \begin{align}
+            \left< \hat{\mathbf{L}} \right>
+            &= \int \phi_a(\mathbf{r}) \left( -i \mathbf{r} \times \nabla \right) \phi_b(\mathbf{r}) d\mathbf{r}\\
+            &= -i
+            \begin{bmatrix}
+            \int \phi_a(\mathbf{r}) y\frac{\partial}{\partial z} \phi_b(\mathbf{r}) d\mathbf{r}
+            - \int \phi_a(\mathbf{r}) z\frac{\partial}{\partial y} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
+            \int \phi_a(\mathbf{r}) z\frac{\partial}{\partial x} \phi_b(\mathbf{r}) d\mathbf{r}
+            - \int \phi_a(\mathbf{r}) x\frac{\partial}{\partial z} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
+            \int \phi_a(\mathbf{r}) x\frac{\partial}{\partial y} \phi_b(\mathbf{r}) d\mathbf{r}
+            - \int \phi_a(\mathbf{r}) y\frac{\partial}{\partial x} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
+            \end{bmatrix}
+        \end{align}
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/electron_repulsion.py b/gbasis/integrals/electron_repulsion.py
index 20b57461..dc3cde8e 100644
--- a/gbasis/integrals/electron_repulsion.py
+++ b/gbasis/integrals/electron_repulsion.py
@@ -208,6 +208,24 @@ def electron_repulsion_integral(
 ):
     """Return the electron repulsion integrals fo the given basis set.
 
+    In the Chemists' notation, the integrals are:
+
+    .. math::
+        
+        \int \phi^*_a(\mathbf{r}_1) \phi_b(\mathbf{r}_1)
+        \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|}
+        \phi^*_c(\mathbf{r}_2) \phi_d(\mathbf{r}_2) d\mathbf{r}
+        
+    And in the Physicists' notation:
+
+    .. math::
+        \int \phi^*_a(\mathbf{r}_1) \phi^*_b(\mathbf{r}_2)
+        \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|}
+        \phi_c(\mathbf{r}_1) \phi_d(\mathbf{r}_2) d\mathbf{r}
+
+    where :math:`\mathbf{r}_1` and :math:`\mathbf{r}_2` are the coordinates of the first and 
+    second electron.
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/kinetic_energy.py b/gbasis/integrals/kinetic_energy.py
index ac67c263..39f7897d 100644
--- a/gbasis/integrals/kinetic_energy.py
+++ b/gbasis/integrals/kinetic_energy.py
@@ -122,6 +122,18 @@ def construct_array_contraction(contractions_one, contractions_two):
 def kinetic_energy_integral(basis, transform=None):
     """Return kinetic energy integral of the given basis set.
 
+    .. math::
+        \begin{align}
+        \left< \hat{T} \right>
+        &= \int \phi_a(\mathbf{r}) \left( -\frac{1}{2} \nabla^2 \right) \phi_b(\mathbf{r}) d\mathbf{r}\\
+        &= -\frac{1}{2}
+        \left(
+        \int \phi_a(\mathbf{r}) \frac{\partial^2}{\partial x^2} \phi_b(\mathbf{r}) d\mathbf{r}
+        + \int \phi_a(\mathbf{r}) \frac{\partial^2}{\partial y^2} \phi_b(\mathbf{r}) d\mathbf{r}
+        + \int \phi_a(\mathbf{r}) \frac{\partial^2}{\partial z^2} \phi_b(\mathbf{r}) d\mathbf{r}
+        \right)
+        \end{align}
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/moment.py b/gbasis/integrals/moment.py
index 44ec36ef..ae9a35de 100644
--- a/gbasis/integrals/moment.py
+++ b/gbasis/integrals/moment.py
@@ -160,6 +160,13 @@ class if the first four indices correspond to the segmented contraction and the
 def moment_integral(basis, moment_coord, moment_orders, transform=None):
     """Return moment integral of the given basis set.
 
+    .. math::
+
+        \int \phi_a (\mathbf{r}) (x - X_C)^{c_x} (y - Y_C)^{c_y} (z - Z_C)^{c_z} \phi_b (\mathbf{r}) d\mathbf{r}
+
+    where :math:`X_C`, :math:`Y_C`, and :math:`Z_C` are the coordinates of the center of the moment, and 
+    :math:`c_x`, :math:`c_y`, and :math:`c_z` are the orders of the moment.
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/momentum.py b/gbasis/integrals/momentum.py
index ac4b1dd9..8e291636 100644
--- a/gbasis/integrals/momentum.py
+++ b/gbasis/integrals/momentum.py
@@ -115,6 +115,18 @@ def construct_array_contraction(contractions_one, contractions_two):
 def momentum_integral(basis, transform=None):
     """Return integral over momentum operator of the given basis set.
 
+    .. math::
+        \begin{align}
+            \left< \hat{\mathbf{p}} \right>
+            &= \int \phi_a(\mathbf{r}) \left( -i \nabla \right) \phi_b(\mathbf{r}) d\mathbf{r}\\
+            &= -i
+            \begin{bmatrix}
+            \int \phi_a(\mathbf{r}) \frac{\partial}{\partial x} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
+            \int \phi_a(\mathbf{r}) \frac{\partial}{\partial y} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
+            \int \phi_a(\mathbf{r}) \frac{\partial}{\partial z} \phi_b(\mathbf{r}) d\mathbf{r}
+            \end{bmatrix}
+        \end{align}
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/nuclear_electron_attraction.py b/gbasis/integrals/nuclear_electron_attraction.py
index 355c2328..a273d31f 100644
--- a/gbasis/integrals/nuclear_electron_attraction.py
+++ b/gbasis/integrals/nuclear_electron_attraction.py
@@ -8,6 +8,12 @@ def nuclear_electron_attraction_integral(
 ):
     """Return the nuclear electron attraction integrals of the basis set in the Cartesian form.
 
+    .. math::
+
+        \int \phi_a(\mathbf{r}) \frac{-Z_C}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
+    
+    where :math:`\mathbf{R}_C` is the coordinate of the nucleus :math:`C` and :math:`Z_C` its charge.
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/overlap.py b/gbasis/integrals/overlap.py
index 5079aaec..f51a6d8b 100644
--- a/gbasis/integrals/overlap.py
+++ b/gbasis/integrals/overlap.py
@@ -111,6 +111,12 @@ def construct_array_contraction(contractions_one, contractions_two):
 def overlap_integral(basis, transform=None):
     """Return overlap integral of the given basis set.
 
+    .. math::
+    
+        \int \phi_a (\mathbf{r}) \phi_b (\mathbf{r}) d\mathbf{r}
+
+    where :math:`\phi_a(\mathbf{r})` and :math:`\phi_b(\mathbf{r})` are the basis functions
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/overlap_asymm.py b/gbasis/integrals/overlap_asymm.py
index bccb23b3..66a2482b 100644
--- a/gbasis/integrals/overlap_asymm.py
+++ b/gbasis/integrals/overlap_asymm.py
@@ -71,6 +71,12 @@ def overlap_integral_asymmetric(
 ):
     """Return overlap integrals between two basis sets.
 
+    .. math::
+    
+        \int \phi_a (\mathbf{r}) \psi_b (\mathbf{r}) d\mathbf{r}
+    
+    where :math:`\phi_a` and :math:`\psi_b` are the basis functions in `basis_one` and `basis_two`.
+
     Parameters
     ----------
     basis_one : list/tuple of GeneralizedContractionShell
diff --git a/gbasis/integrals/point_charge.py b/gbasis/integrals/point_charge.py
index 7a615696..a71783ac 100644
--- a/gbasis/integrals/point_charge.py
+++ b/gbasis/integrals/point_charge.py
@@ -271,6 +271,12 @@ def point_charge_integral(
 ):
     r"""Return the point-charge interaction integrals of basis set in the given coordinate systems.
 
+    .. math::
+
+        \int \phi_a(\mathbf{r}) \frac{1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
+    
+    where :math:`\mathbf{R}_C` is the position of the point charge.
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell