diff --git a/doc/source/code/figure/closed_system.py.diff b/doc/source/code/figure/closed_system.py.diff
index 6ad93be562a7e793ca1f825fa8af2ec5641ef114..748e9eca48693b42d4f54faf64e988e034096101 100644
--- a/doc/source/code/figure/closed_system.py.diff
+++ b/doc/source/code/figure/closed_system.py.diff
@@ -60,14 +60,9 @@
  #HIDDEN_BEGIN_yvri
  def plot_spectrum(syst, Bfields):
  
-     # In the following, we compute the spectrum of the quantum dot
-     # using dense matrix methods. This works in this toy example, as
-     # the system is tiny. In a real example, one would want to use
-     # sparse matrix methods
- 
      energies = []
      for B in Bfields:
-         # Obtain the Hamiltonian as a dense matrix
+         # Obtain the Hamiltonian as a sparse matrix
          ham_mat = syst.hamiltonian_submatrix(params=dict(B=B), sparse=True)
  
          # we only calculate the 15 lowest eigenvalues
diff --git a/doc/source/pre/whatsnew/1.4.rst b/doc/source/pre/whatsnew/1.4.rst
index aba6608d1f4ea3aa9b0b8eaf759c704edf5e4d50..566c92df07bc550500c368eb619ba9e0416b34ce 100644
--- a/doc/source/pre/whatsnew/1.4.rst
+++ b/doc/source/pre/whatsnew/1.4.rst
@@ -83,7 +83,9 @@ Here is an example for extracting the symmetry group of a graphene system::
 
 ``kwant.qsymm`` also contains functionality for converting Qsymm models to Kwant Builders,
 and vice versa, and for working with continuum Hamiltonians (such as would be used with
-``kwant.continuum``)
+``kwant.continuum``).
+This integration requires separately installing Qsymm, which is available on the
+`Python Package Index <https://pypi.org/project/qsymm/>`_.
 
 Automatic Peierls phase calculation
 -----------------------------------
@@ -101,8 +103,7 @@ which calculates the Peierls phases for you::
       return -t * peierls(a, b)
 
   syst = make_system(hopping)
-  lead = make_lead(hopping)
-  lead.substituted(peierls='peierls_lead')
+  lead = make_lead(hopping).substituted(peierls='peierls_lead')
   syst.attach_lead(lead)
   syst = syst.finalized()
 
@@ -135,19 +136,19 @@ has a parameter ``V``, and one wishes to have different values for ``V`` in
 the scattering region and leads, one could do the following::
 
    syst = kwant.Builder()
-   syst.fill(model.substituted(V='V_dot', ...))
+   syst.fill(model.substituted(V='V_dot'), ...))
 
    lead = kwant.Builder()
    lead.fill(model.substituted(V='V_lead'), ...)
 
    syst.attach_lead(lead)
-   fsyst = syst.finalized()
+   syst = syst.finalized()
 
    kwant.smatrix(syst, params=dict(V_dot=0, V_lead=1))
 
 Interpolated density plots
 --------------------------
-A new function, `~kwant.plotter.density`, has been added that can be used to
+A new function, `kwant.plotter.density`, has been added that can be used to
 visualize a density defined over the sites of a Kwant system. This convolves
 the "discrete" density (defined over the system sites) with a "bump" function
 in realspace. The output of `~kwant.plotter.density` can be more informative
@@ -166,7 +167,8 @@ but in an inconsistent way (e.g. a parameter 'phi' that is a superconducting
 phase in one value function, but a peierls phase in another). This leads
 to bugs that are confusing and hard to track down.
 
-Concretely, the above means that the following no longer works::
+For this reason value functions may no longer have default values for paramters.
+Concretely this means that the following no longer works::
 
   syst = kwant.Builder()
 
@@ -190,7 +192,7 @@ To deal with many parameters, the following idiom may be useful::
   ...
   smatrix = kwant.smatrix(syst, E, params=dict(defaults, d=4, e=5))
 
-Note that it allows to override defaults as well as to add additional
+Note that this allows to override defaults as well as to add additional
 parameters.
 
 System parameters can now be inspected
@@ -228,9 +230,9 @@ This is a provisional API that may be changed in a future version of Kwant.
 Passing system arguments via ``args`` is deprecated in favor of ``params``
 --------------------------------------------------------------------------
 It is now deprecated to pass arguments to systems by providing the
-``args`` parameter (in ``kwant.smatrix`` and elsewhere). This is
-error prone and requires that all value functions take the same
-formal parameters, even if they do not depend on all of them. The
+``args`` parameter (in ``kwant.smatrix`` and elsewhere). Passing arguments
+via ``args`` is error prone and requires that all value functions take the
+same formal parameters, even if they do not depend on all of them. The
 preferred way of passing parameters to Kwant systems is by passing
 a dictionary using ``params``::
 
@@ -244,7 +246,7 @@ a dictionary using ``params``::
   # Compare this to the deprecated 'args'
   kwant.smatrix(syst, args=(0.5, 0.2))
 
-The ability to provide ``args`` will be removed in a future Kwant version.
+Providing ``args`` will be removed in a future Kwant version.
 
 Finalized Builders keep track of which sites were added when attaching leads
 ----------------------------------------------------------------------------
@@ -261,8 +263,8 @@ that this option is not available in `~kwant.plotter.current`.  In order to use
 it, one has to call ``streamplot`` directly as shown in the docstring of
 ``current``.
 
-kwant.continuum.discretize can be used with rectangular lattices
-----------------------------------------------------------------
+`kwant.continuum.discretize` can be used with rectangular lattices
+------------------------------------------------------------------
 Previously the discretizer could only be used with lattices with the same
 lattice constant in all directions. Now it is possible to pass rectangular
 lattices to the discretizer::
diff --git a/doc/source/reference/kwant.kpm.rst b/doc/source/reference/kwant.kpm.rst
index adec65d795005c7e1aaab2ac2931446f49cfeef6..877bf03d042e5e2c90141caaba8e542b5ce2d8a5 100644
--- a/doc/source/reference/kwant.kpm.rst
+++ b/doc/source/reference/kwant.kpm.rst
@@ -4,4 +4,4 @@
 .. automodule:: kwant.kpm
    :members:
    :special-members:
-   :exclude-members: __weakref__
+   :exclude-members: __weakref__, __init__
diff --git a/doc/source/reference/kwant.physics.rst b/doc/source/reference/kwant.physics.rst
index 5c764d7dd796298d84203bca28fa0cfe06137d29..a2d361c7ac81390e71d1344d3e9908df09a390e0 100644
--- a/doc/source/reference/kwant.physics.rst
+++ b/doc/source/reference/kwant.physics.rst
@@ -27,5 +27,6 @@ Computation of magnetic field gauge
 
 .. autosummary::
    :toctree: generated/
+   :template: autosummary/functor.rst
 
    magnetic_gauge
diff --git a/doc/source/tutorial/discretize.rst b/doc/source/tutorial/discretize.rst
index b5edb0157f8a4ba5c3542b525328a48d852a8d42..b737b00ca074cbfd88a7e863be67ebc03063d902 100644
--- a/doc/source/tutorial/discretize.rst
+++ b/doc/source/tutorial/discretize.rst
@@ -161,8 +161,8 @@ energy eigenstates:
 .. image:: /code/figure/discretizer_gs.*
 
 Note in the above that we pass the spatially varying potential *function*
-to our system via a parameter called ``V``, because the symbol $V$
-was used in the intial, symbolic, definition of the Hamiltonian.
+to our system via a parameter called ``V``, because the symbol :math:`V`
+was used in the initial, symbolic, definition of the Hamiltonian.
 
 In addition, the function passed as ``V`` expects two input parameters ``x``
 and ``y``, the same as in the initial continuum Hamiltonian.
diff --git a/doc/source/tutorial/faq.rst b/doc/source/tutorial/faq.rst
index 9c1871e7675e757a3d868e48adda607638e364cc..c2ba0700dbc8d3ef4aaf3c2665fbf4ed77a944eb 100644
--- a/doc/source/tutorial/faq.rst
+++ b/doc/source/tutorial/faq.rst
@@ -423,7 +423,7 @@ the modes are sorted in the following way:
     + Positive velocity modes are ordered by *decreasing* momentum
 
 For more complicated systems and band structures this can lead to some
-possibly unintuitive orderings:
+unintuitive orderings:
 
 .. image:: /code/figure/faq_pm2.*
 
@@ -437,9 +437,9 @@ infinity" (*not* towards the system) which means that the incoming modes are
 those that have *negative* velocities.
 
 This means that for a lead attached on the left of a scattering region (with
-symmetry vector $(-1, 0)$, for example), the
-positive $k$ direction (when inspecting the lead's band structure) actually
-corresponds to the *negative* $x$ direction.
+symmetry vector :math:`(-1, 0)`, for example), the
+positive :math:`k` direction (when inspecting the lead's band structure) actually
+corresponds to the *negative* :math:`x` direction.
 
 
 How does Kwant order components of an individual wavefunction?
diff --git a/doc/source/tutorial/graphene.rst b/doc/source/tutorial/graphene.rst
index b4281a80f370128c5c9860384c33ae6a20a9cb85..63bcf5f29859ea09f11b84410c64871bf49fe9b5 100644
--- a/doc/source/tutorial/graphene.rst
+++ b/doc/source/tutorial/graphene.rst
@@ -112,9 +112,6 @@ in the following piece of code:
     :start-after: #HIDDEN_BEGIN_zydk
     :end-before: #HIDDEN_END_zydk
 
-Here we use in contrast to the previous example a sparse matrix and
-the sparse linear algebra functionality of SciPy.
-
 The code for computing the band structure and the conductance is identical
 to the previous examples, and needs not be further explained here.
 
diff --git a/doc/source/tutorial/kpm.rst b/doc/source/tutorial/kpm.rst
index b8e7bd2378f57bf1c10cb5bd5803e41c936e0b41..38c84891cbd84ea044e8600ee56c3e3d0a1cd638 100644
--- a/doc/source/tutorial/kpm.rst
+++ b/doc/source/tutorial/kpm.rst
@@ -340,7 +340,7 @@ vectors
 Note that the Kubo conductivity must be normalized with the area covered
 by the vectors. In this case, each local vector represents a site, and
 covers an area of half a unit cell, while the sum covers one unit cell.
-It is possible to use random vectors to get an average spectation
+It is possible to use random vectors to get an average expectation
 value of the conductivity over large parts of the system. In this
 case, the area that normalizes the result, is the area covered by
 the random vectors.
diff --git a/doc/source/tutorial/spectrum.rst b/doc/source/tutorial/spectrum.rst
index 46843440019be37e7bc0d3e99124366eb6df68ee..734f5a1d0e161d2434789bc9ce3f3a880e122419 100644
--- a/doc/source/tutorial/spectrum.rst
+++ b/doc/source/tutorial/spectrum.rst
@@ -122,7 +122,7 @@ better for the special case of a square lattice.
 
 
 As our model breaks time reversal symmetry (because of the applied magnetic
-field) we can also see an intereting property of the eigenstates, namely
+field) we can also see an interesting property of the eigenstates, namely
 that they can have *non-zero* local current. We can calculate the local
 current due to a state by using `kwant.operator.Current` and plotting
 it using `kwant.plotter.current`:
diff --git a/kwant/kpm.py b/kwant/kpm.py
index 8cf8db4d063ae7d72254382b094423c864fa5dd1..967ba4f144fa6bde127d71c9d89ff76255790c0d 100644
--- a/kwant/kpm.py
+++ b/kwant/kpm.py
@@ -41,7 +41,7 @@ class SpectralDensity:
     .. math::
        ρ_A(E) = ρ(E) A(E),
 
-    where :math:`ρ(E) = \\sum_{k=0}^{D-1} δ(E-E_k)` is the density of
+    where :math:`ρ(E) = \sum_{k=0}^{D-1} δ(E-E_k)` is the density of
     states, and :math:`A(E)` is the expectation value of :math:`A` for
     all the eigenstates with energy :math:`E`.
 
@@ -253,11 +253,9 @@ class SpectralDensity:
 
         Returns
         -------
-        ``(energies, densities)`` if the ``energy`` parameter is not
-        provided, else ``densities``.
-
         energies : array of floats
             Drawn from the nodes of the highest Chebyshev polynomial.
+            Not returned if 'energy' was not provided
         densities : float or array of floats
             single ``float`` if the ``energy`` parameter is a single
             ``float``, else an array of ``float``.
@@ -593,7 +591,7 @@ class Correlator:
         self._build_integral_factor()
 
     def __call__(self, mu=0, temperature=0):
-        """Returns the linear response χ_{α β}(µ, T)
+        """Returns the linear response :math:`χ_{α β}(µ, T)`
 
         Parameters
         ----------
diff --git a/kwant/physics/gauge.py b/kwant/physics/gauge.py
index 1fb80cfa1e3538f6f5ca151abb11a38c147ebb09..bb2b3bc6d9949743a5fe915d57fbe75d371c65f0 100644
--- a/kwant/physics/gauge.py
+++ b/kwant/physics/gauge.py
@@ -944,7 +944,7 @@ class magnetic_gauge:
 
     Parameters
     ----------
-    syst : kwant.builder.FiniteSystem or kwant.builder.InfiniteSystem
+    syst : `kwant.builder.FiniteSystem` or `kwant.builder.InfiniteSystem`
 
     Examples
     --------
@@ -1027,7 +1027,7 @@ class magnetic_gauge:
             The first callable computes the Peierls phase in the scattering
             region and the remaining callables compute the Peierls phases
             in each of the leads. Each callable takes a pair of
-            `~kwant.builder.Site`s 'a, b' and returns a unit complex
+            `~kwant.builder.Site` (a hopping) and returns a unit complex
             number (Peierls phase) that multiplies that hopping.
         """
         return self._peierls(syst_field, *lead_fields, tol=tol, average=False)