diff --git a/doc/source/tutorial/discretize.rst b/doc/source/tutorial/discretize.rst
index 31c4b3bb3ca217888c3c0e1b53e444a270906705..9417d20b52e8e6448f32e64894928f3eefa56abc 100644
--- a/doc/source/tutorial/discretize.rst
+++ b/doc/source/tutorial/discretize.rst
@@ -202,7 +202,7 @@ energy eigenstates:
 
         ham = syst.hamiltonian_submatrix(params=dict(V=potential), sparse=True)
         evecs = scipy.sparse.linalg.eigsh(ham, k=10, which='SM')[1]
-        kwant.plotter.map(syst, abs(evecs[:, n])**2, show=False)
+        kwant.plotter.density(syst, abs(evecs[:, n])**2, show=False)
 
 .. jupyter-execute::
     :hide-code:
@@ -281,7 +281,7 @@ and plot its dispersion using `kwant.plotter.bands`:
     pyplot.show()
 
 In the above we see the edge states of the quantum spin Hall effect, which
-we can visualize using `kwant.plotter.map`:
+we can visualize using `kwant.plotter.density`:
 
 .. jupyter-execute::
 
@@ -298,8 +298,8 @@ we can visualize using `kwant.plotter.map`:
     rho_sz = sum(spin_density(psi) for psi in wf(0))  # states from left lead
 
     fig, (ax1, ax2) = pyplot.subplots(1, 2, sharey=True, figsize=(16, 4))
-    kwant.plotter.map(syst, wf_sqr, ax=ax1)
-    kwant.plotter.map(syst, rho_sz, ax=ax2)
+    kwant.plotter.density(syst, wf_sqr, ax=ax1)
+    kwant.plotter.density(syst, rho_sz, ax=ax2)
 
     ax = ax1
     im = [obj for obj in ax.get_children()
diff --git a/doc/source/tutorial/kpm.rst b/doc/source/tutorial/kpm.rst
index 7e1232581335aa1a77969899f59344e34c93eba7..f1d6650c44a2ae2b08016e74ffa357e6b643caf4 100644
--- a/doc/source/tutorial/kpm.rst
+++ b/doc/source/tutorial/kpm.rst
@@ -196,7 +196,7 @@ object that represents the density of states for this system.
 
     fsyst = make_syst().finalized()
 
-    spectrum = kwant.kpm.SpectralDensity(fsyst)
+    spectrum = kwant.kpm.SpectralDensity(fsyst, rng=0)
 
 The `~kwant.kpm.SpectralDensity` can then be called like a function to obtain a
 sequence of energies in the spectrum of the Hamiltonian, and the corresponding
@@ -319,7 +319,8 @@ and plot their respective local density of states.
     # 'num_vectors' can be unspecified when using 'LocalVectors'
     local_dos = kwant.kpm.SpectralDensity(fsyst_staggered, num_vectors=None,
                                           vector_factory=vector_factory,
-                                          mean=False)
+                                          mean=False,
+                                          rng=0)
     energies, densities = local_dos()
 
 .. jupyter-execute::
@@ -346,7 +347,7 @@ The simplest way to obtain a more accurate solution is to use the
 .. jupyter-execute::
     :hide-code:
 
-    spectrum = kwant.kpm.SpectralDensity(fsyst)
+    spectrum = kwant.kpm.SpectralDensity(fsyst, rng=0)
     original_dos = spectrum()
 
 .. jupyter-execute::
@@ -399,7 +400,7 @@ shape as the system Hamiltonian.
 
     # identity matrix
     matrix_op = scipy.sparse.eye(len(fsyst.sites))
-    matrix_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=matrix_op)
+    matrix_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=matrix_op, rng=0)
 
 Or, to do the same calculation using `kwant.operator.Density`:
 
@@ -407,7 +408,7 @@ Or, to do the same calculation using `kwant.operator.Density`:
 
     # 'sum=True' means we sum over all the sites
     kwant_op = kwant.operator.Density(fsyst, sum=True)
-    operator_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op)
+    operator_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op, rng=0)
 
 Spectral density with random vectors
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -420,7 +421,7 @@ sum over all the sites of the system:
 
     # 'sum=False' is the default, but we include it explicitly here for clarity.
     kwant_op = kwant.operator.Density(fsyst, sum=False)
-    local_dos = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op)
+    local_dos = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op, rng=0)
 
 `~kwant.kpm.SpectralDensity` will properly handle this vector output,
 and will average the local density obtained with random vectors.
@@ -499,11 +500,13 @@ vectors
     # component 'xx'
     s_factory = kwant.kpm.LocalVectors(fsyst_topo, where)
     cond_xx = kwant.kpm.conductivity(fsyst_topo, alpha='x', beta='x', mean=True,
-                                     num_vectors=None, vector_factory=s_factory)
+                                     num_vectors=None, vector_factory=s_factory,
+                                     rng=0)
     # component 'xy'
     s_factory = kwant.kpm.LocalVectors(fsyst_topo, where)
     cond_xy = kwant.kpm.conductivity(fsyst_topo, alpha='x', beta='y', mean=True,
-                                     num_vectors=None, vector_factory=s_factory)
+                                     num_vectors=None, vector_factory=s_factory,
+                                     rng=0)
 
     energies = cond_xx.energies
     cond_array_xx = np.array([cond_xx(e, temperature=0.01) for e in energies])
@@ -528,7 +531,8 @@ the random vectors.
 
     s_factory = kwant.kpm.LocalVectors(fsyst_topo, where)
     spectrum = kwant.kpm.SpectralDensity(fsyst_topo, num_vectors=None,
-                                         vector_factory=s_factory)
+                                         vector_factory=s_factory,
+                                         rng=0)
 
     plot_dos_and_curves(
     (spectrum.energies, spectrum.densities * 8),
@@ -570,7 +574,8 @@ and which returns a vector in that Hilbert space:
            yield np.rint(np.random.random_sample(n)) * 2 - 1
 
     custom_factory = kwant.kpm.SpectralDensity(fsyst,
-                                               vector_factory=binary_vectors())
+                                               vector_factory=binary_vectors(),
+                                               rng=0)
 
 Aditionally, a `~kwant.kpm.LocalVectors` generator is also available, that
 returns local vectors that correspond to the sites passed. Note that
@@ -613,8 +618,8 @@ methods for computing the local density of states, one using
     def rho_alt(bra, ket):
         return np.vdot(bra, ket)
 
-    rho_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho)
-    rho_alt_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho_alt)
+    rho_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho, rng=0)
+    rho_alt_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho_alt, rng=0)
 
 __ operator_spectral_density_
 
diff --git a/doc/source/tutorial/operators.rst b/doc/source/tutorial/operators.rst
index 1646ae86e2cd5a47747843f6dd952144b0a92175..897c919c4b0b77e82b5ad39efec28c664df3d75a 100644
--- a/doc/source/tutorial/operators.rst
+++ b/doc/source/tutorial/operators.rst
@@ -165,7 +165,7 @@ inside the scattering region. The z component is shown by the color scale:
     def plot_densities(syst, densities):
         fig, axes = plt.subplots(1, len(densities), figsize=(13, 10))
         for ax, (title, rho) in zip(axes, densities):
-            kwant.plotter.map(syst, rho, ax=ax, a=4)
+            kwant.plotter.density(syst, rho, ax=ax)
             ax.set_title(title)
         plt.show()