diff --git a/doc/source/images/2-spin_orbit.py.diff b/doc/source/images/2-spin_orbit.py.diff
index d935aa9411b4f4d7c86871e931c7c6228d3508f0..7e151abe363e2f2bf54e111348288d0ec475a8bb 100644
--- a/doc/source/images/2-spin_orbit.py.diff
+++ b/doc/source/images/2-spin_orbit.py.diff
@@ -3,11 +3,11 @@
 @@ -17,6 +17,7 @@
  
  # For matrix support
- import numpy
+ import tinyarray
 +import latex, html
  
  # define Pauli-matrices for convenience
- sigma_0 = numpy.eye(2)
+ sigma_0 = tinyarray.array([[1, 0], [0, 1]])
 @@ -73,19 +74,24 @@
          smatrix = kwant.solve(sys, energy)
          data.append(smatrix.transmission(1, 0))
diff --git a/doc/source/images/5-superconductor_band_structure.py.diff b/doc/source/images/5-superconductor_band_structure.py.diff
index 7392cd639bca0e8f253f26978dedfb8a30aaea65..d99af4979c1b62c592d487f207bcf0955e9c2161 100644
--- a/doc/source/images/5-superconductor_band_structure.py.diff
+++ b/doc/source/images/5-superconductor_band_structure.py.diff
@@ -1,14 +1,14 @@
 --- original
 +++ modified
-@@ -16,6 +16,7 @@
+@@ -17,6 +17,7 @@
  
  # For plotting
  from matplotlib import pyplot
 +import latex, html
  
- tau_x = np.array([[0, 1], [1, 0]])
- tau_z = np.array([[1, 0], [0, -1]])
-@@ -46,12 +47,19 @@
+ tau_x = tinyarray.array([[0, 1], [1, 0]])
+ tau_z = tinyarray.array([[1, 0], [0, -1]])
+@@ -47,12 +48,19 @@
      # the bandstructure
      energy_list = [lead.energies(k) for k in momenta]
  
diff --git a/doc/source/tutorial/2-spin_orbit.py b/doc/source/tutorial/2-spin_orbit.py
index 18c225935f272abe6fe24e1eabdab70b320ab684..d52996f365b52733236575757ce13c56b0a3d55a 100644
--- a/doc/source/tutorial/2-spin_orbit.py
+++ b/doc/source/tutorial/2-spin_orbit.py
@@ -17,15 +17,15 @@ from matplotlib import pyplot
 
 # For matrix support
 #HIDDEN_BEGIN_xumz
-import numpy
+import tinyarray
 #HIDDEN_END_xumz
 
 # define Pauli-matrices for convenience
 #HIDDEN_BEGIN_hwbt
-sigma_0 = numpy.eye(2)
-sigma_x = numpy.array([[0, 1], [1, 0]])
-sigma_y = numpy.array([[0, -1j], [1j, 0]])
-sigma_z = numpy.array([[1, 0], [0, -1]])
+sigma_0 = tinyarray.array([[1, 0], [0, 1]])
+sigma_x = tinyarray.array([[0, 1], [1, 0]])
+sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
+sigma_z = tinyarray.array([[1, 0], [0, -1]])
 #HIDDEN_END_hwbt
 
 
diff --git a/doc/source/tutorial/5-superconductor_band_structure.py b/doc/source/tutorial/5-superconductor_band_structure.py
index db1c99850a16583cfaf545fce83b90c51c164976..e26c6407e807566bbaa1e2353e37fcd1c1a38e14 100644
--- a/doc/source/tutorial/5-superconductor_band_structure.py
+++ b/doc/source/tutorial/5-superconductor_band_structure.py
@@ -13,13 +13,14 @@
 import kwant
 
 import numpy as np
+import tinyarray
 
 # For plotting
 from matplotlib import pyplot
 
 #HIDDEN_BEGIN_nbvn
-tau_x = np.array([[0, 1], [1, 0]])
-tau_z = np.array([[1, 0], [0, -1]])
+tau_x = tinyarray.array([[0, 1], [1, 0]])
+tau_z = tinyarray.array([[1, 0], [0, -1]])
 
 
 def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
diff --git a/doc/source/tutorial/tutorial2.rst b/doc/source/tutorial/tutorial2.rst
index f78e0d576efcb90f6b80357f7bbe3ca33517199c..4bb79ac7b7ad8a617597cf20cc367bd82b44fb65 100644
--- a/doc/source/tutorial/tutorial2.rst
+++ b/doc/source/tutorial/tutorial2.rst
@@ -31,11 +31,11 @@ we will show that a very simple extension of our previous examples will
 exactly show this behavior (Note though that no care was taken to choose
 realistic parameters).
 
-The tight-binding model corresponding to the Rashba-Hamiltonian
-naturally exhibits a 2x2-matrix structure of onsite energies and hoppings.
-In order to deal with matrices in python, kwant uses the `numpy package
-<numpy.scipy.org>`_. In order to use matrices in our program, we thus also
-have to import that package:
+The tight-binding model corresponding to the Rashba-Hamiltonian naturally
+exhibits a 2x2-matrix structure of onsite energies and hoppings.  In order to
+use matrices in our program, we import the tinyarray package.  (`NumPy
+<http://numpy.scipy.org/>`_ would work as well, but tinyarray is much faster
+for small arrays.)
 
 .. literalinclude:: 2-spin_orbit.py
     :start-after: #HIDDEN_BEGIN_xumz
@@ -96,6 +96,22 @@ the following, clearly non-monotonic conductance steps:
 
 .. specialnote:: Technical details
 
+  - The tinyarray package, one of the dependencies of kwant, implements
+    efficient small arrays.  It is used internally in kwant for storing small
+    vectors and matrices.  For performance, it is preferable to define small
+    arrays that are going to be used with kwant using tinyarray.  However,
+    NumPy would work as well::
+
+        import numpy
+        sigma_0 = numpy.array([[1, 0], [0, 1]])
+        sigma_x = numpy.array([[0, 1], [1, 0]])
+        sigma_y = numpy.array([[0, -1j], [1j, 0]])
+        sigma_z = numpy.array([[1, 0], [0, -1]])
+
+    tinyarray arrays behave for most purposes like NumPy arrays except that
+    they are immutable: they cannot be changed once created.  This is important
+    for kwant: it allows them to be used directly as dictionary keys.
+
   - It should be emphasized that the relative hopping used for
     `~kwant.builder.Builder.possible_hoppings` is given in terms of
     lattice indices, i.e. relative to the Bravais lattice vectors.