diff --git a/Docs/source/diffusion.rst b/Docs/source/diffusion.rst
index bf5bf1165b..3fb847da68 100644
--- a/Docs/source/diffusion.rst
+++ b/Docs/source/diffusion.rst
@@ -4,7 +4,7 @@
 Thermal Diffusion
 *****************
 
-Castro incorporates explicit thermal diffusion into the energy equation.
+Castro incorporates explicit thermal diffusion into the energy equations.
 In terms of the specific internal energy, :math:`e`, this appears as:
 
 .. math:: \rho \frac{De}{Dt} + p \nabla \cdot \ub = \nabla \cdot \kth \nabla T
@@ -20,15 +20,82 @@ where :math:`\kth` is the thermal conductivity, with units
 
      USE_DIFFUSION=TRUE
 
-It is treated explicitly, by constructing the contribution to the evolution as a
-source term. This is time-centered to achieve second-order accuracy
+Thermal Diffusion related source codes are contained in the ``diffusion`` directory.
+Thermal Diffusion is treated explicitly, by constructing the contribution to the
+evolution as a source term. This is time-centered to achieve second-order accuracy
 in time.
 
+Overall Procedure
+=================
+
+Computing Thermal Conductivity
+------------------------------
+The main function that computes the diffusion term is ``getTempDiffusionTerm()``.
+Within ``getTempDiffusionTerm()``, it first calculates the cell centered
+thermal conductivity, :math:`\kth` contained in the variable ``coeff_cc``
+using the function ``fill_temp_cond()`` located in ``diffusion_util.cpp``.
+``fill_temp_cond()`` fills an ``eos_state`` using the
+input conserved variables, which is used to calculate :math:`\kth` via
+``conductivity(eos_state)``. ``conductivity()`` routine is supplied via
+the ``Microphysics`` package. See :ref:`sec:conductivities` to see the
+specific choices of conductivity routines available.
+
+.. note::
+   The diffusion approximation breaks down at the surface of stars,
+   where the density rapidly drops and the mean free path becomes
+   large. In those instances, you should use the flux limited diffusion
+   module in Castro to evolve a radiation field.
+
+Now :math:`\kth` is reset to 0 unless
+:math:`\rho \gt \mathrm{castro::diffuse\_cutoff\_density}`.
+And if :math:`\rho \lt \mathrm{castro::diffuse\_cutoff\_density\_hi}`,
+a linear scaling of :math:`\kth` is done as:
+
+.. math::
+
+   \kth = \kth \cdot \frac{\rho - \mathtt{castro.diffuse\_cutoff\_density}}{\mathtt{castro.diffuse\_cutoff\_density\_hi} - \mathtt{castro.diffuse\_cutoff\_density}}
+
+Lastly, :math:`\kth` is scaled with ``castro::diffuse_cond_scale_fac``,
+a runtime parameter controlled by the user.
+
+After obtaining cell-centered :math:`\kth`, we do an average along
+i, j, and k depending on the direction to obtain face-centered MultiFabs.
+This is stored in ``coeffs``, a vector of MultiFabs, and the number of
+MultiFabs corresponds to geometry dimension, since a :math:`\nabla` operator
+will be applied to it later.
+These Multifabs have 1 ghost cells due to the nature of MLMG solvers.
+
+.. _sec:thermal_diffusion:
+
+Computing Thermal Diffusion
+---------------------------
+We are now ready to compute :math:`\nabla \cdot \kth \nabla T`
+after obtaining :math:`\kth`. This is done in the ``applyop_mlmg()`` function
+in ``Diffusion.cpp``. It defines ``mlabec`` an instance of class
+``MLABecLaplacian`` which defines the Laplacian of the form:
+
+.. math::
+   (A\alpha - B\nabla \cdot \beta \nabla) \phi = f
+
+where A and B are constant scalars, and :math:`\alpha` and :math:`\beta`
+are scalar fields. In order to make it correctly represents our diffusion term,
+we set A = 0 and B = -1, which is done via ``mlabec.setScalars(0.0, -1.0)``.
+Now we recognize :math:`\beta = \kth`, which needs to be an array of MultiFab,
+corresponding to dimension. This is done via ``mlabec.setBCoeffs()``.
+
+One of the important flag that we need to pass in is to set ``setMetricTerm(true)``.
+This enables modifications due to curvilinear coordinates.
+
+Finally we create an instance of ``MLMG`` using ``mlabec``, and call
+``mlmg.apply()``, which simply evaluates the LHS but do not solve it.
+See more information in the amrex documentation:
+https://amrex-codes.github.io/amrex/docs_html/LinearSolvers.html
+
 
 Timestep Limiter
 ================
 
-Castro integrates diffusion explicitly in time—this means that
+Castro integrates diffusion explicitly in time; this means that
 there is a diffusion timestep limiter.
 
 To see the similarity to the thermal diffusion equation, consider the
@@ -59,27 +126,16 @@ The following parameter affects diffusion:
 
       castro.diffuse_temp = 1
       castro.do_hydro = 0
+* ``castro.diffuse_cond_scale_fac``: a linear scaling to :math:`\kth`. (default 0).
 
-.. index:: castro.diffusion_cutoff_density, castro.diffusion_cutoff_density_hi
-
-The diffusion approximation breaks down at the surface of stars,
-where the density rapidly drops and the mean free path becomes
-large. In those instances, you should use the flux limited diffusion
-module in Castro to evolve a radiation field. However, if your
-interest is only on the diffusion in the interior, you can use
-the parameters:
+* ``castro.diffuse_cutoff_density``: density under which :math:`\kth` is set to 0.
+  (Default: -1e200)
 
- * ``castro.diffuse_cutoff_density``
+* ``castro.diffuse_cutoff_density_hi``: density under which a linear scaling is
+  applied to :math:`\kth`, see section :ref:`sec:thermal_diffusion` for details.
+  (Default: -1e200)
 
- * ``castro.diffuse_cutoff_density_hi``
-
-to specify a density,
-below which, diffusion is not modeled. This is implemented in the
-code by linearly scaling the conductivity to zero between these limits, e.g.,
-
-.. math::
-
-   \kth = \kth \cdot \frac{\rho - \mathtt{castro.diffuse\_cutoff\_density}}{\mathtt{castro.diffuse\_cutoff\_density\_hi} - \mathtt{castro.diffuse\_cutoff\_density}}
+.. _sec:conductivities:
 
 Conductivities
 ==============