diff --git a/Docs/sphinx_documentation/source/LinearSolvers.rst b/Docs/sphinx_documentation/source/LinearSolvers.rst
index 6a08c0dc5a7..b617c686faf 100644
--- a/Docs/sphinx_documentation/source/LinearSolvers.rst
+++ b/Docs/sphinx_documentation/source/LinearSolvers.rst
@@ -35,11 +35,12 @@ below
                      const Vector<BoxArray>& a_grids,
                      const Vector<DistributionMapping>& a_dmap,
                      const LPInfo& a_info = LPInfo(),
-                     const Vector<FabFactory<FArrayBox> const*>& a_factory = {});
+                     const Vector<FabFactory<FArrayBox> const*>& a_factory = {},
+                     const int a_ncomp = 1);
 
 It takes :cpp:`Vectors` of :cpp:`Geometry`, :cpp:`BoxArray` and
 :cpp:`DistributionMapping`.  The arguments are :cpp:`Vectors` because MLMG can
-do multi-level composite solve.  If you are using it for single-level,
+do multi-level composite solves.  If you are using it for single-level,
 you can do
 
 .. highlight:: c++
@@ -60,25 +61,72 @@ After the linear operator is built, we need to set up boundary
 conditions.  This will be discussed later in section
 :ref:`sec:linearsolver:bc`.
 
-For :cpp:`MLABecLaplacian`, we next need to call member functions
+Coefficients
+------------
+
+Next, we consider the coefficients for equation :eq:`eqn::abeclap`.
+For :cpp:`MLPoisson`, there are no coefficients to set so nothing needs to be done.
+For :cpp:`MLABecLaplacian`, we need to call member functions :cpp:`setScalars`,
+:cpp:`setACoeffs`, and :cpp:`setBCoeffs`.
+The :cpp:`setScalars` function sets the scalar constants :math:`A` and :math:`B`
+
+.. code-block::
+
+    void setScalars (Real a, Real b) noexcept;
+
+
+For the general case where
+:math:`\alpha` and :math:`\beta` are scalar fields, we use
+
+  .. code-block::
+
+      void setACoeffs (int amrlev, const MultiFab& alpha);
+      void setBCoeffs (int amrlev, const Array<MultiFab const*,AMREX_SPACEDIM>& beta);
+
+For the case where :math:`\alpha` and/or :math:`\beta` are scalar constants,
+there is the option to use
+
+  .. code-block::
+
+      void setACoeffs (int amrlev, Real alpha);
+      void setBCoeffs (int amrlev, Real beta);
+      void setBCoeffs (int amrlev, Vector<Real> const& beta);
+
+Note, however, that the solver behaviour is the same regardless of which functions you
+use to set the coefficients. These functions solely copy the constant value(s) to a MultiFab
+internal to ``MLMG`` and so no appreciable efficiency gains can be expected.
+
+For :cpp:`MLNodeLaplacian`,
+one can set a variable :cpp:`sigma` with the member function
 
 .. highlight:: c++
 
 ::
 
-    void setScalars (Real A, Real B);
-    void setACoeffs (int amrlev, const MultiFab& alpha);
-    void setBCoeffs (int amrlev, const Array<MultiFab const*,AMREX_SPACEDIM>& beta);
+    void setSigma (int amrlev, const MultiFab& a_sigma);
 
-to set up the coefficients for equation :eq:`eqn::abeclap`. This is unnecessary for
-:cpp:`MLPoisson`, as there are no coefficients to set.  For :cpp:`MLNodeLaplacian`,
-one needs to call the member function
+or a constant :cpp:`sigma` during declaration or definition
 
 .. highlight:: c++
 
 ::
 
-    void setSigma (int amrlev, const MultiFab& a_sigma);
+    MLNodeLaplacian (const Vector<Geometry>& a_geom,
+                     const Vector<BoxArray>& a_grids,
+                     const Vector<DistributionMapping>& a_dmap,
+                     const LPInfo& a_info = LPInfo(),
+                     const Vector<FabFactory<FArrayBox> const*>& a_factory = {},
+                     Real  a_const_sigma = Real(0.0));
+
+    void define (const Vector<Geometry>& a_geom,
+                 const Vector<BoxArray>& a_grids,
+                 const Vector<DistributionMapping>& a_dmap,
+                 const LPInfo& a_info = LPInfo(),
+                 const Vector<FabFactory<FArrayBox> const*>& a_factory = {},
+                 Real  a_const_sigma = Real(0.0));
+
+Here, setting a constant :cpp:`sigma` alters the internal behavior of the solver making it more
+efficient for this special case.
 
 The :cpp:`int amrlev` parameter should be zero for single-level
 solves.  For multi-level solves, each level needs to be provided with