diff --git a/src/openfermion/chem/reduced_hamiltonian.py b/src/openfermion/chem/reduced_hamiltonian.py
index 40dedb06c..98a58d1ef 100644
--- a/src/openfermion/chem/reduced_hamiltonian.py
+++ b/src/openfermion/chem/reduced_hamiltonian.py
@@ -15,10 +15,11 @@ def make_reduced_hamiltonian(
     lift the 1-body terms to the two-body space.
 
     Derivation:
-        use the fact that i^l = (1/(n -1)) sum_{jk}\delta_{jk}i^ j^ k l
-                          i^l = (-1/(n -1)) sum_{jk}\delta_{jk}j^ i^ k l
-                          i^l = (-1/(n -1)) sum_{jk}\delta_{jk}i^ j^ l k
-                          i^l = (1/(n -1)) sum_{jk}\delta_{jk}j^ i^ l k
+        use the fact that
+        - $i^l = (1/(n -1)) sum_{jk}\delta_{jk}i^ j^ k l$
+        - $i^l = (-1/(n -1)) sum_{jk}\delta_{jk}j^ i^ k l$
+        - $i^l = (-1/(n -1)) sum_{jk}\delta_{jk}i^ j^ l k$
+        - $i^l = (1/(n -1)) sum_{jk}\delta_{jk}j^ i^ l k$
 
         Rewrite each one-body term as an even weighting of all four 2-RDM
         elements with delta functions. Then rearrange terms so that each ijkl
diff --git a/src/openfermion/circuits/gates/fermionic_simulation.py b/src/openfermion/circuits/gates/fermionic_simulation.py
index fd8213bde..3c651717b 100644
--- a/src/openfermion/circuits/gates/fermionic_simulation.py
+++ b/src/openfermion/circuits/gates/fermionic_simulation.py
@@ -52,6 +52,7 @@ def state_swap_eigen_component(x: str, y: str, sign: int = 1, angle: float = 0):
     """The +/- eigen-component of the operation that swaps states x and y.
 
     For example, state_swap_eigen_component('01', '10', ±1) with angle θ returns
+
         ┌                               ┐
         │0, 0,           0,            0│
         │0, 0.5,         ±0.5 e^{-iθ}, 0│
@@ -102,7 +103,7 @@ def fermionic_simulation_gates_from_interaction_operator(
 ):
     r"""
     Given $H = \sum_{I \subset [n]} H_I$, returns gates
-    $\left\{G_I\right\} = \left\{e^{i H_I\right\}$.
+    $\left\{G_I\right\} = \left\{e^{i H_I}\right\}$.
 
     Each term $H_I$ is the sum of all terms in $H$ that involve exactly the
     orbitals $I$.
diff --git a/src/openfermion/circuits/primitives/bogoliubov_transform.py b/src/openfermion/circuits/primitives/bogoliubov_transform.py
index 7b7fcb985..be6e6d8b8 100644
--- a/src/openfermion/circuits/primitives/bogoliubov_transform.py
+++ b/src/openfermion/circuits/primitives/bogoliubov_transform.py
@@ -51,7 +51,7 @@ def bogoliubov_transform(
     given by
 
     $$
-        b^\dagger_p = \sum_{q=1}^N W_{pq} a^\dagger_q
+        b^\dagger_p = \sum_{q=1}^N W_{pq} a^\dagger_q \\
                       + \sum_{q=N+1}^{2N} W_{pq} a_q.
     $$
 
diff --git a/src/openfermion/hamiltonians/__init__.py b/src/openfermion/hamiltonians/__init__.py
index 5ad591c44..0de50d665 100644
--- a/src/openfermion/hamiltonians/__init__.py
+++ b/src/openfermion/hamiltonians/__init__.py
@@ -56,6 +56,7 @@
 
 from .special_operators import (
     s_plus_operator,
+    s_minus_operator,
     s_squared_operator,
     sx_operator,
     sy_operator,
diff --git a/src/openfermion/hamiltonians/general_hubbard.py b/src/openfermion/hamiltonians/general_hubbard.py
index 37d195231..50d2a2a85 100644
--- a/src/openfermion/hamiltonians/general_hubbard.py
+++ b/src/openfermion/hamiltonians/general_hubbard.py
@@ -63,6 +63,7 @@ class FermiHubbardModel:
     form
 
     $$
+        \begin{align}
         H = &- \sum_{a < b} t_{a, b}^{(\mathrm{onsite})}
                \sum_{i} \sum_{\sigma}
                      (a^\dagger_{i, a, \sigma} a_{i, b, \sigma} +
@@ -104,6 +105,7 @@ class FermiHubbardModel:
             \\
             &- h \sum_{i} \sum_{a}
                 \left(n_{i, a, \uparrow} - n_{i, a, \downarrow}\right)
+        \end{align}
     $$
 
     where
@@ -159,6 +161,7 @@ class FermiHubbardModel:
             &- \sum_{a} \mu_a
                \sum_i n_{i, a}
         \end{align}
+    $$
     """
 
     def __init__(
diff --git a/src/openfermion/hamiltonians/richardson_gaudin.py b/src/openfermion/hamiltonians/richardson_gaudin.py
index 4de5e6a72..7c99aa9c2 100644
--- a/src/openfermion/hamiltonians/richardson_gaudin.py
+++ b/src/openfermion/hamiltonians/richardson_gaudin.py
@@ -22,30 +22,30 @@ class RichardsonGaudin(DOCIHamiltonian):
     r"""Richardson Gaudin model.
 
     Class for storing and constructing Richardson Gaudin hamiltonians
-    combining an equi-distant potential ladder like potential per
+    combining an equidistant potential ladder like potential per
     qubit with a uniform coupling between any pair of
-    qubits with coupling strength g, which can be either attractive
-    (g<0) or repulsive (g>0).
+    qubits with coupling strength `g`, which can be either attractive
+    (`g<0`) or repulsive (`g>0`).
 
-    The operators represented by this class has the form:
-
-    .. math::
+    The operators represented by this class have the form:
 
+    $$
         H = \sum_{p=0} (p + 1) N_p + g/2 \sum_{p < q} P_p^\dagger P_q,
+    $$
 
     where
 
-    .. math::
-
+    $$
         \begin{align}
         N_p &= (1 - \sigma^Z_p)/2, \\
         P_p &= a_{p,\beta} a_{p,\alpha} = S^{-} = \sigma^X + i \sigma^Y, \\
         g &= constant coupling term
         \end{align}
+    $$
 
     Note;
         The diagonal of the Hamiltonian is composed of the values in
-        range((n_qubits+1)*n_qubits//2+1).
+        `range((n_qubits+1)*n_qubits//2+1)`.
     """
 
     def __init__(self, g, n_qubits):
diff --git a/src/openfermion/hamiltonians/special_operators.py b/src/openfermion/hamiltonians/special_operators.py
index d4b4ffd72..fa235a55d 100644
--- a/src/openfermion/hamiltonians/special_operators.py
+++ b/src/openfermion/hamiltonians/special_operators.py
@@ -48,7 +48,7 @@ def s_plus_operator(n_spatial_orbitals: int) -> FermionOperator:
 
 
 def s_minus_operator(n_spatial_orbitals: int) -> FermionOperator:
-    r"""Return the s+ operator.
+    r"""Return the s- operator.
 
     $$
         \begin{align}
diff --git a/src/openfermion/linalg/erpa.py b/src/openfermion/linalg/erpa.py
index 33f89f83a..b93cbb057 100644
--- a/src/openfermion/linalg/erpa.py
+++ b/src/openfermion/linalg/erpa.py
@@ -12,7 +12,7 @@ def erpa_eom_hamiltonian(
     h_ijkl: numpy.ndarray, tpdm: numpy.ndarray, p: int, q: int, r: int, s: int
 ) -> Union[float, complex]:
     """
-    Evaluate sum_{a,b,c,d}h_{a, b, d, c}<psi[p^ q, [a^ b^ c d, r^ s]]psi>
+    Evaluate $\sum_{a,b,c,d}h_{a, b, d, c}<\psi[p^ q, [a^ b^ c d, r^ s]]\psi>$
 
     Args:
         h_ijkl: two-body integral tensors of full reduced Hamiltonian
@@ -84,7 +84,7 @@ def singlet_erpa(
     """
     Generate the singlet ERPA equations
 
-    [ea + eb, [H, sa, sb]] = [ea, [H, sa]]
+        [ea + eb, [H, sa, sb]] = [ea, [H, sa]]
 
     The erpa equations are solved with scipy.linalg.eig which calls
     lapack's geev
diff --git a/src/openfermion/linalg/rdm_reconstruction.py b/src/openfermion/linalg/rdm_reconstruction.py
index da3a2a456..6bd01cc04 100644
--- a/src/openfermion/linalg/rdm_reconstruction.py
+++ b/src/openfermion/linalg/rdm_reconstruction.py
@@ -6,7 +6,7 @@ def valdemoro_reconstruction(tpdm, n_electrons):
     """
     Build a 3-RDM by cumulant expansion and setting 3rd cumulant to zero
 
-    d3 approx = D ^ D ^ D + 3 (2C) ^ D
+        d3 approx = D ^ D ^ D + 3 (2C) ^ D
 
     tpdm has normalization (n choose 2) where n is the number of electrons
 
diff --git a/src/openfermion/linalg/sparse_tools.py b/src/openfermion/linalg/sparse_tools.py
index 453db13d3..fcc93257d 100644
--- a/src/openfermion/linalg/sparse_tools.py
+++ b/src/openfermion/linalg/sparse_tools.py
@@ -49,8 +49,9 @@ def jordan_wigner_ladder_sparse(n_qubits, tensor_factor, ladder_type):
     r"""Make a matrix representation of a fermion ladder operator.
 
     Operators are mapped as follows:
-    a_j^\dagger -> Z_0 .. Z_{j-1} (X_j - iY_j) / 2
-    a_j -> Z_0 .. Z_{j-1} (X_j + iY_j) / 2
+
+        a_j^\dagger -> Z_0 .. Z_{j-1} (X_j - iY_j) / 2
+        a_j -> Z_0 .. Z_{j-1} (X_j + iY_j) / 2
 
     Args:
         index: This is a nonzero integer. The integer indicates the tensor
@@ -75,8 +76,9 @@ def jordan_wigner_sparse(fermion_operator, n_qubits=None):
     r"""Initialize a Scipy sparse matrix from a FermionOperator.
 
     Operators are mapped as follows:
-    a_j^\dagger -> Z_0 .. Z_{j-1} (X_j - iY_j) / 2
-    a_j -> Z_0 .. Z_{j-1} (X_j + iY_j) / 2
+
+        a_j^\dagger -> Z_0 .. Z_{j-1} (X_j - iY_j) / 2
+        a_j -> Z_0 .. Z_{j-1} (X_j + iY_j) / 2
 
     Args:
         fermion_operator(FermionOperator): instance of the FermionOperator
diff --git a/src/openfermion/linalg/wave_fitting.py b/src/openfermion/linalg/wave_fitting.py
index 2f186647c..4cff2ad06 100644
--- a/src/openfermion/linalg/wave_fitting.py
+++ b/src/openfermion/linalg/wave_fitting.py
@@ -43,7 +43,7 @@ def prony(signal: numpy.ndarray) -> Tuple[numpy.ndarray, numpy.ndarray]:
     """Estimates amplitudes and phases of a sparse signal using Prony's method.
 
     Single-ancilla quantum phase estimation returns a signal
-    g(k)=sum (aj*exp(i*k*phij)), where aj and phij are the amplitudes
+    $g(k)=\sum (aj*exp(i*k*phij))$, where aj and phij are the amplitudes
     and corresponding eigenvalues of the unitary whose phases we wish
     to estimate. When more than one amplitude is involved, Prony's method
     provides a simple estimation tool, which achieves near-Heisenberg-limited
diff --git a/src/openfermion/linalg/wedge_product.py b/src/openfermion/linalg/wedge_product.py
index ab4d2f55b..154332a13 100644
--- a/src/openfermion/linalg/wedge_product.py
+++ b/src/openfermion/linalg/wedge_product.py
@@ -73,14 +73,14 @@ def wedge(left_tensor, right_tensor, left_index_ranks, right_index_ranks):
     The wedge product is defined as
 
     $$
-        \\begin{align}
-        a_{j_{1}, j_{2}, ...,j_{p}}^{i_{1}, i_{2}, ..., i_{p}} \\wedge
+        \begin{align}
+        a_{j_{1}, j_{2}, ...,j_{p}}^{i_{1}, i_{2}, ..., i_{p}} \wedge
         b_{j_{p+1}, j_{p+2}, ..., j_{N}}^{i_{p+1}, i_{p + 2}, ..., i_{N}} =
-        \\left(\\frac{1}{N!}\\right)^{2} = \\sum_{\\pi, \\sigma}\\epsilon(\\pi)
-        \\epsilon(\\sigma)\\pi \\sigma
+        \left(\frac{1}{N!}\right)^{2} = \sum_{\pi, \sigma}\epsilon(\pi)
+        \epsilon(\sigma)\pi \sigma
         a_{j_{1}, j_{2}, ...,j_{p}}^{i_{1}, i_{2}, ..., i_{p}}
         b_{j_{p+1}, j_{p+2}, ..., j_{N}}^{i_{p+1}, i_{p + 2}, ..., i_{N}}
-        \\end{align}
+        \end{align}
     $$
 
     The top indices are those that transform contravariently.  The bottom
diff --git a/src/openfermion/measurements/fermion_partitioning.py b/src/openfermion/measurements/fermion_partitioning.py
index 62e514bff..e4d67b9c0 100644
--- a/src/openfermion/measurements/fermion_partitioning.py
+++ b/src/openfermion/measurements/fermion_partitioning.py
@@ -22,16 +22,16 @@ def pair_within(labels: list) -> list:
 
     A pairing of a list is a set of pairs of list elements. E.g. a pairing of
 
-    labels = [1, 2, 3, 4, 5, 6, 7, 8]
+        labels = [1, 2, 3, 4, 5, 6, 7, 8]
 
     could be
 
-    [(1, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 2), (3, 4), (5, 6), (7, 8)]
 
     (Note that we insist each element only appears in a pairing once; the
     following is not a pairing:
 
-    [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
 
     This function generates a set of len(labels)-1 such pairings of the input
     list labels, such that each element in labels is paired with each other
@@ -86,16 +86,16 @@ def pair_between(frag1: list, frag2: list, start_offset: int = 0) -> tuple:
 
     A pairing of a list is a set of pairs of list elements. E.g. a pairing of
 
-    labels = [1, 2, 3, 4, 5, 6, 7, 8]
+        labels = [1, 2, 3, 4, 5, 6, 7, 8]
 
     could be
 
-    [(1, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 2), (3, 4), (5, 6), (7, 8)]
 
     (Note that we insist each element only appears in a pairing once; the
     following is not a pairing:
 
-    [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
 
     This function generates a set of pairings between elements of frag1
     and frag2 such that element1 in frag 1 and element2 in frag2,
@@ -202,16 +202,16 @@ def pair_within_simultaneously(labels: list) -> tuple:
 
     A pairing of a list is a set of pairs of list elements. E.g. a pairing of
 
-    labels = [1, 2, 3, 4, 5, 6, 7, 8]
+        labels = [1, 2, 3, 4, 5, 6, 7, 8]
 
     could be
 
-    [(1, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 2), (3, 4), (5, 6), (7, 8)]
 
     (Note that we insist each element only appears in a pairing once; the
     following is not a pairing:
 
-    [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
 
     This function generates a set of pairings such that for every four elements
     (i,j,k,l) in 'labels', there exists one pairing containing both (i,j) and
@@ -372,16 +372,16 @@ def pair_within_simultaneously_binned(binned_majoranas: list) -> tuple:
 
     A pairing of a list is a set of pairs of list elements. E.g. a pairing of
 
-    labels = [1, 2, 3, 4, 5, 6, 7, 8]
+        labels = [1, 2, 3, 4, 5, 6, 7, 8]
 
     could be
 
-    [(1, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 2), (3, 4), (5, 6), (7, 8)]
 
     (Note that we insist each element only appears in a pairing once; the
     following is not a pairing:
 
-    [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
 
     This function generates a pairing of a list of Majoranas that covers
     all 2-RDM elements that conserve a set of symmetry conditions. That is,
@@ -431,16 +431,16 @@ def pair_within_simultaneously_symmetric(num_fermions: int, num_symmetries: int)
 
     A pairing of a list is a set of pairs of list elements. E.g. a pairing of
 
-    labels = [1, 2, 3, 4, 5, 6, 7, 8]
+        labels = [1, 2, 3, 4, 5, 6, 7, 8]
 
     could be
 
-    [(1, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 2), (3, 4), (5, 6), (7, 8)]
 
     (Note that we insist each element only appears in a pairing once; the
     following is not a pairing:
 
-    [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
+        [(1, 1), (2, 2), (3, 4), (5, 6), (7, 8)]
 
     This function generates a pairing of a list of Majoranas that covers
     all 2-RDM elements that conserve a set of symmetry conditions. That is,
diff --git a/src/openfermion/measurements/rdm_equality_constraints.py b/src/openfermion/measurements/rdm_equality_constraints.py
index ec3facdfa..edfc6da69 100644
--- a/src/openfermion/measurements/rdm_equality_constraints.py
+++ b/src/openfermion/measurements/rdm_equality_constraints.py
@@ -16,10 +16,10 @@
 def one_body_fermion_constraints(n_orbitals, n_fermions):
     """Generates one-body positivity constraints on fermionic RDMs.
 
-        The specific constraints implemented are known positivity constraints
-        on the one-fermion reduced density matrices. Constraints are generated
-        in the form of FermionOperators whose expectation value is known to be
-        zero for any N-Representable state. Generators are used for efficiency.
+    The specific constraints implemented are known positivity constraints
+    on the one-fermion reduced density matrices. Constraints are generated
+    in the form of FermionOperators whose expectation value is known to be
+    zero for any N-Representable state. Generators are used for efficiency.
 
     Args:
         n_orbitals(int): number of spin-orbitals on which operators act.
@@ -48,10 +48,10 @@ def one_body_fermion_constraints(n_orbitals, n_fermions):
 def two_body_fermion_constraints(n_orbitals, n_fermions):
     """Generates two-body positivity constraints on fermionic RDMs.
 
-        The specific constraints implemented are known positivity constraints
-        on the two-fermion reduced density matrices. Constraints are generated
-        in the form of FermionOperators whose expectation value is known to be
-        zero for any N-Representable state. Generators are used for efficiency.
+    The specific constraints implemented are known positivity constraints
+    on the two-fermion reduced density matrices. Constraints are generated
+    in the form of FermionOperators whose expectation value is known to be
+    zero for any N-Representable state. Generators are used for efficiency.
 
     Args:
         n_orbitals(int): number of spin-orbitals on which operators act.
diff --git a/src/openfermion/ops/operators/quad_operator.py b/src/openfermion/ops/operators/quad_operator.py
index 98eeada4b..16f294d51 100644
--- a/src/openfermion/ops/operators/quad_operator.py
+++ b/src/openfermion/ops/operators/quad_operator.py
@@ -18,8 +18,10 @@ class QuadOperator(SymbolicOperator):
     """QuadOperator stores a sum of products of canonical quadrature operators.
 
     They are defined in terms of the bosonic ladder operators:
-    q = sqrt{hbar/2}(b+b^)
-    p = -isqrt{hbar/2}(b-b^)
+
+        q = sqrt{hbar/2}(b+b^)
+        p = -isqrt{hbar/2}(b-b^)
+
     where hbar is a constant appearing in the commutator of q and p:
     [q, p] = i hbar
 
diff --git a/src/openfermion/ops/representations/interaction_operator.py b/src/openfermion/ops/representations/interaction_operator.py
index 32a384eeb..4287c3e84 100644
--- a/src/openfermion/ops/representations/interaction_operator.py
+++ b/src/openfermion/ops/representations/interaction_operator.py
@@ -34,7 +34,7 @@ class InteractionOperator(PolynomialTensor):
     class take the form:
 
         $$
-            constant + \sum_{p, q} h_{p, q} a^\dagger_p a_q +
+            \text{constant} + \sum_{p, q} h_{p, q} a^\dagger_p a_q +
             \sum_{p, q, r, s} h_{p, q, r, s} a^\dagger_p a^\dagger_q a_r a_s.
         $$
 
diff --git a/src/openfermion/resource_estimates/pbc/thc/thc_integrals.py b/src/openfermion/resource_estimates/pbc/thc/thc_integrals.py
index da323fc62..c262cfe30 100644
--- a/src/openfermion/resource_estimates/pbc/thc/thc_integrals.py
+++ b/src/openfermion/resource_estimates/pbc/thc/thc_integrals.py
@@ -70,14 +70,14 @@ def __init__(
     def get_eri(self, ikpts: list) -> npt.NDArray:
         r"""Construct ERIs given kpt indices.
 
-        .. math::
-
+        $$
             (pkp qkq| rkr sks) = \\sum_{mu nu} zeta[iq, dG, dG', mu, nu]
             chi[kp,p,mu]* chi[kq,q,mu] chi[kp,p,nu]* chi[ks,s,nu]
+        $$
 
         Arguments:
             ikpts: list of four integers representing the index of the kpoint in
-                self.kmf.kpts
+                `self.kmf.kpts`
 
         Returns:
             eris: ([pkp][qkq]|[rkr][sks])
diff --git a/src/openfermion/transforms/opconversions/bksf.py b/src/openfermion/transforms/opconversions/bksf.py
index 00c090baf..35d4ecce6 100644
--- a/src/openfermion/transforms/opconversions/bksf.py
+++ b/src/openfermion/transforms/opconversions/bksf.py
@@ -53,11 +53,19 @@ def bravyi_kitaev_fast_interaction_op(iop: InteractionOperator) -> QubitOperator
     The electronic Hamiltonian is represented in terms of creation and
     annihilation operators. These creation and annihilation operators could be
     used to define Majorana modes as follows:
+
+    $$
         c_{2i} = a_i + a^{\dagger}_i,
         c_{2i+1} = (a_i - a^{\dagger}_{i})/(1j)
+    $$
+
     These Majorana modes can be used to define edge operators B_i and A_{ij}:
+
+    $$
         B_i=c_{2i}c_{2i+1},
         A_{ij}=c_{2i}c_{2j}
+    $$
+
     using these edge operators the fermionic algebra can be generated and
     hence all the terms in the electronic Hamiltonian can be expressed in
     terms of edge operators. The terms in electronic Hamiltonian can be
diff --git a/src/openfermion/transforms/opconversions/conversions.py b/src/openfermion/transforms/opconversions/conversions.py
index 95f45b63c..9dbf49ebe 100644
--- a/src/openfermion/transforms/opconversions/conversions.py
+++ b/src/openfermion/transforms/opconversions/conversions.py
@@ -166,6 +166,7 @@ def get_majorana_operator(
 
     Uses the convention of even + odd indexing of Majorana modes derived from
     a fermionic mode:
+
         fermion annhil.  c_k  -> ( gamma_{2k} + 1.j * gamma_{2k+1} ) / 2
         fermion creation c^_k -> ( gamma_{2k} - 1.j * gamma_{2k+1} ) / 2
 
diff --git a/src/openfermion/transforms/opconversions/jordan_wigner.py b/src/openfermion/transforms/opconversions/jordan_wigner.py
index 60131e709..17f649b48 100644
--- a/src/openfermion/transforms/opconversions/jordan_wigner.py
+++ b/src/openfermion/transforms/opconversions/jordan_wigner.py
@@ -25,8 +25,11 @@ def jordan_wigner(operator):
     to a QubitOperator.
 
     Operators are mapped as follows:
-    a_j^\dagger -> Z_0 .. Z_{j-1} (X_j - iY_j) / 2
-    a_j -> Z_0 .. Z_{j-1} (X_j + iY_j) / 2
+
+    $$
+        a_j^\dagger -> Z_0 .. Z_{j-1} (X_j - iY_j) / 2
+        a_j -> Z_0 .. Z_{j-1} (X_j + iY_j) / 2
+    $$
 
     Returns:
         transformed_operator: An instance of the QubitOperator class.
diff --git a/src/openfermion/transforms/opconversions/reverse_jordan_wigner.py b/src/openfermion/transforms/opconversions/reverse_jordan_wigner.py
index fbd32012c..e340c5edf 100644
--- a/src/openfermion/transforms/opconversions/reverse_jordan_wigner.py
+++ b/src/openfermion/transforms/opconversions/reverse_jordan_wigner.py
@@ -21,9 +21,12 @@ def reverse_jordan_wigner(qubit_operator, n_qubits=None):
     Jordan-Wigner transform.
 
     Operators are mapped as follows:
-    Z_j -> I - 2 a^\dagger_j a_j
-    X_j -> (a^\dagger_j + a_j) Z_{j-1} Z_{j-2} .. Z_0
-    Y_j -> i (a^\dagger_j - a_j) Z_{j-1} Z_{j-2} .. Z_0
+
+    $$
+        Z_j -> I - 2 a^\dagger_j a_j
+        X_j -> (a^\dagger_j + a_j) Z_{j-1} Z_{j-2} .. Z_0
+        Y_j -> i (a^\dagger_j - a_j) Z_{j-1} Z_{j-2} .. Z_0
+    $$
 
     Args:
         qubit_operator: the QubitOperator to be transformed.
diff --git a/src/openfermion/transforms/repconversions/qubit_operator_transforms.py b/src/openfermion/transforms/repconversions/qubit_operator_transforms.py
index d5161abce..ff83ff052 100644
--- a/src/openfermion/transforms/repconversions/qubit_operator_transforms.py
+++ b/src/openfermion/transforms/repconversions/qubit_operator_transforms.py
@@ -120,7 +120,7 @@ def rotate_qubit_by_pauli(qop, pauli, angle):
     r"""
     Rotate qubit operator by exponential of Pauli.
 
-    Perform the rotation e^{-i \theta * P}Qe^{i \theta * P}
+    Perform the rotation $e^{-i \theta * P}Qe^{i \theta * P}$
     on a qubitoperator Q and a Pauli operator P.
 
     Args:
diff --git a/src/openfermion/transforms/repconversions/weyl_ordering.py b/src/openfermion/transforms/repconversions/weyl_ordering.py
index 1df907d3c..33dada554 100644
--- a/src/openfermion/transforms/repconversions/weyl_ordering.py
+++ b/src/openfermion/transforms/repconversions/weyl_ordering.py
@@ -44,8 +44,8 @@ def weyl_polynomial_quantization(polynomial):
     The Weyl quantization is performed by applying McCoy's formula
     directly to a polynomial term of the form q^m p^n:
 
-    q^m p^n ->
-        (1/ 2^n) sum_{r=0}^{n} Binomial(n, r) \hat{q}^r \hat{p}^m q^{n-r}
+    $q^m p^n ->
+        (1/ 2^n) sum_{r=0}^{n} Binomial(n, r) \hat{q}^r \hat{p}^m q^{n-r}$
 
     where q and p are phase space variables, and \hat{q} and \hat{p}
     are quadrature operators.
@@ -111,7 +111,7 @@ def symmetric_ordering(operator, ignore_coeff=True, ignore_identity=True):
     The symmetric ordering is performed by applying McCoy's formula
     directly to polynomial terms of quadrature operators:
 
-    q^m p^n -> (1/ 2^n) sum_{r=0}^{n} Binomial(n, r) q^r p^m q^{n-r}
+    $q^m p^n -> (1/ 2^n) sum_{r=0}^{n} Binomial(n, r) q^r p^m q^{n-r}$
 
     Note: in general, symmetric ordering is performed on a single term
     containing the tensor product of various operators. However, this
diff --git a/src/openfermion/utils/bch_expansion.py b/src/openfermion/utils/bch_expansion.py
index 8e75748ea..d2a2f1770 100644
--- a/src/openfermion/utils/bch_expansion.py
+++ b/src/openfermion/utils/bch_expansion.py
@@ -16,7 +16,7 @@
 
 
 def bch_expand(*ops, **kwargs):
-    """Compute log[e^{x_1} ... e^{x_N}] using the BCH formula.
+    """Compute $\log[e^{x_1} ... e^{x_N}]$ using the BCH formula.
 
     This implementation is explained in arXiv:1712.01348.
 
@@ -65,7 +65,7 @@ def _bch_expand_multiple_terms(*ops, **kwargs):
 
 
 def _bch_expand_two_terms(x, y, order=6):
-    """Compute log[e^x e^y] using the Baker-Campbell-Hausdorff formula.
+    """Compute $\log[e^x e^y]$ using the Baker-Campbell-Hausdorff formula.
 
     Args:
         x: An operator for which multiplication and addition are supported.