diff --git a/doc/source/images/FAQ.py.diff b/doc/source/images/FAQ.py.diff
index b8649809d26ecbb1279ce6b781ce2bef4c13e45c..4608a56f3651b2be0aa62a20b941421e2cc9d24e 100644
--- a/doc/source/images/FAQ.py.diff
+++ b/doc/source/images/FAQ.py.diff
@@ -1,223 +1,204 @@
 --- original
 +++ modified
-@@ -1,4 +1,4 @@
--
+@@ -1,6 +1,7 @@
+ # Frequently Asked Questions
+ # ==========================
+ 
 +import _defs
- from cmath import exp
- import numpy as np
  import kwant
-@@ -11,6 +11,13 @@
+ import numpy as np
+ import tinyarray
+@@ -9,6 +10,12 @@
+ matplotlib.rcParams['figure.figsize'] = (3.5, 3.5)
  
- ######## What is a Site? ##############
  
 +def save_figure(file_name):
-+    if not file_name:
-+        return
 +    for extension in ('pdf', 'png'):
-+        plt.savefig('.'.join((file_name,extension)),
++        plt.savefig('.'.join((file_name, extension)),
 +                    dpi=_defs.dpi, bbox_inches='tight')
 +
++
+ ######## What is a Site? ##############
+ 
  a = 1
- lat = kwant.lattice.square(a)
- syst = kwant.Builder()
 @@ -19,6 +26,7 @@
- syst[lat(1,1)] = 4
+ syst[lat(1, 1)] = 4
  
  kwant.plot(syst)
-+save_figure('FAQ122')
++save_figure("FAQ122")
  
  
  
-@@ -40,7 +48,14 @@
- syst[lat2(0,0)] = 4 ## syst[subB(0,0)] = 4
+@@ -40,6 +48,7 @@
+ syst[lat2(0, 0)] = 4  # syst[subB(0, 0)] = 4
  
  kwant.plot(syst)
-+save_figure('FAQ123')
- 
-+def save_figure(file_name):
-+    if not file_name:
-+        return
-+    for extension in ('pdf', 'png'):
-+        plt.savefig('.'.join((file_name,extension)),
-+                    dpi=_defs.dpi, bbox_inches='tight')
++save_figure("FAQ123")
  
  
  ########## What is a hopping? #######################
-@@ -54,6 +69,7 @@
- syst[(lat(1,0) , lat(1,1))] = -1
+@@ -53,6 +62,7 @@
+ syst[(lat(1, 0), lat(1, 1))] = -1
  
  kwant.plot(syst)
-+save_figure('FAQ124')
- 
++save_figure("FAQ124")
  
+ ########### How to make a hole in a system? ###################"
  
-@@ -73,7 +89,7 @@
+@@ -69,6 +79,7 @@
  
- syst[ (lat(i,j,k) for i in range(L) for j in range(W) for k in range(H)) ] = 4
+ syst[(lat(i, j, k) for i in range(L) for j in range(W) for k in range(H))] = 4
  kwant.plot(syst)
--
-+save_figure('FAQ1')
++save_figure("FAQ1")
  
- ## Deleting sites to create a hole
+ # Delete sites to create a hole
  
-@@ -84,6 +100,7 @@
-             if ((L-2)*a/4 <= x <= 3*L*a/4) and (W*a/4 <= y <= 3*(W-2)*a/4) and (0 <= z <= (H-1)*a):
-                 del syst[lat(i,j,k)]
- kwant.plot(syst)
-+save_figure('FAQ2')
+@@ -80,6 +91,7 @@
+     del syst[site]
  
+ kwant.plot(syst)
++save_figure("FAQ2")
  
  
-@@ -145,7 +162,7 @@
+ ################ How can we get access to the sites of our system? ####################
+@@ -133,7 +145,7 @@
  
  ## Plotting the system
  plot_system(syst)
 -
-+save_figure('FAQ6B')
++save_figure("FAQ6B")
  
  
+ ############### How to create every hoppings in a given direction using Hoppingkind? ################
+@@ -149,6 +161,7 @@
  
-@@ -162,7 +179,7 @@
- 
- syst[kwant.builder.HoppingKind((1,0), lat)] = -1
+ syst[kwant.builder.HoppingKind((1, 0), lat)] = -1
  kwant.plot(syst)
--
-+save_figure('FAQ3')
- ## Polyatomic lattice
++save_figure("FAQ3")
  
- lat = kwant.lattice.kagome()
-@@ -193,12 +210,13 @@
+ # Polyatomic lattice
  
- syst[kwant.builder.HoppingKind((0,1), b, b)] = -1 # equivalent to syst[kwant.builder.HoppingKind((0,1), b)] = -1
+@@ -181,12 +194,14 @@
+ # equivalent to syst[kwant.builder.HoppingKind((0, 1), b)] = -1
+ syst[kwant.builder.HoppingKind((0, 1), b, b)] = -1
  plot_system(syst)
-+save_figure('FAQ10')
- del syst[kwant.builder.HoppingKind((0,1), b, b)] ## delete the hoppings previously created
- syst[kwant.builder.HoppingKind((0,0), a, b)] = -1
- syst[kwant.builder.HoppingKind((0,0), a, c)] = -1
- syst[kwant.builder.HoppingKind((0,0), c, b)] = -1
++save_figure("FAQ10")
+ # Delete the hoppings previously created
+ del syst[kwant.builder.HoppingKind((0, 1), b, b)]
+ syst[kwant.builder.HoppingKind((0, 0), a, b)] = -1
+ syst[kwant.builder.HoppingKind((0, 0), a, c)] = -1
+ syst[kwant.builder.HoppingKind((0, 0), c, b)] = -1
  plot_system(syst)
--
-+save_figure('FAQ11')
++save_figure("FAQ11")
  
  
  ########## How to create the hoppings between adjacent sites? ################
-@@ -213,11 +231,13 @@
+@@ -201,11 +216,13 @@
  
- syst[lat.neighbors()] = -1 ## equivalent to lat.neighbors(1)
+ syst[lat.neighbors()] = -1  # Equivalent to lat.neighbors(1)
  kwant.plot(syst)
-+save_figure('FAQ4')
++save_figure("FAQ4")
  
- del syst[lat.neighbors()] ## deletes every hoppings previously created to add new one
+ del syst[lat.neighbors()]  # Delete all nearest-neighbor hoppings
  syst[lat.neighbors(2)] = -1
  
  kwant.plot(syst)
-+save_figure('FAQ5')
++save_figure("FAQ5")
  
- ## Polyatomic lattice
+ # Polyatomic lattice
+ 
+@@ -224,14 +241,17 @@
  
-@@ -236,14 +256,17 @@
- ## .neighbors()
  syst[lat.neighbors()] = -1
  plot_system(syst)
-+save_figure('FAQ7')
- del syst[lat.neighbors()] ## delete the hoppings previously created
++save_figure("FAQ7")
+ del syst[lat.neighbors()]  # Delete the hoppings previously created
  syst[a.neighbors()] = -1
  plot_system(syst)
-+save_figure('FAQ8')
- del syst[a.neighbors()] ## deletes every hoppings previously created to add new one
++save_figure("FAQ8")
+ del syst[a.neighbors()]  # Delete the hoppings previously created
  
  
  syst[lat.neighbors(2)] = -1
  plot_system(syst)
-+save_figure('FAQ9')
++save_figure("FAQ9")
  del syst[lat.neighbors(2)]
  
  
-@@ -282,9 +305,10 @@
- syst[ (subB(i,j) for i in range(L) for j in range(W)) ] = 4
+@@ -264,6 +284,7 @@
+ syst[(subB(i, j) for i in range(L) for j in range(W))] = 4
  syst[lat.neighbors()] = -1
  plot_system(syst)
-+save_figure('FAQAB')
- ## We manually add sites of the same lead lattice
- 
--
-+a = 2
- lat2 = kwant.lattice.square(a)
- 
- def shapetop(pos):
-@@ -303,6 +327,7 @@
- syst[((lat2(i,-1), subA(i,0)) for i in range(5))] = -1
- syst[((lat2(i+2,5), subB(i,4)) for i in range(5))] = -1
++save_figure("FAQAA")
+ 
+ # Create a lead
+ lat_lead = kwant.lattice.square()
+@@ -273,6 +294,7 @@
+ lead1[(lat_lead(i, 0) for i in range(2, 7))] = 4
+ lead1[lat_lead.neighbors()] = -1
+ plot_system(lead1)
++save_figure("FAQAB")
+ 
+ syst[(lat_lead(i, 5) for i in range(2, 7))] = 4
+ syst[lat_lead.neighbors()] = -1
+@@ -280,9 +302,11 @@
+ # Manually attach sites from graphene to square lattice
+ syst[((lat_lead(i+2, 5), subB(i, 4)) for i in range(5))] = -1
  plot_system(syst)
-+save_figure('FAQAC')
- ## Creation of the top lead
- 
- lat_lead = kwant.lattice.square(a)
-@@ -320,6 +345,8 @@
++save_figure("FAQAC")
  
  syst.attach_lead(lead1)
  plot_system(syst)
-+save_figure('FAQAD')
-+
- ## Creation of the bottom lead
- 
- sym_lead2 = kwant.TranslationalSymmetry((0,-a))
-@@ -335,6 +362,7 @@
- syst.attach_lead(lead2)
- 
- plot_system(syst)
-+save_figure('FAQAE')
- 
++save_figure("FAQAD")
  
  
-@@ -356,6 +384,7 @@
+ ############# How to cut a finite system out of a system with translationnal symmetries? ###########
+@@ -302,6 +326,7 @@
+ cuboid = kwant.Builder()
  cuboid.fill(model, cuboid_shape, (0, 0, 0));
- 
  kwant.plot(cuboid);
-+save_figure('FAQaaa')
++save_figure("FAQaaa")
  
  # Build electrode (black).
  def electrode_shape(site):
-@@ -369,7 +398,7 @@
+@@ -316,6 +341,7 @@
  
  cuboid.attach_lead(electrode)
  kwant.plot(cuboid);
--
-+save_figure('FAQbbb')
++save_figure("FAQbbb")
  
- ###### How to extract the wavefunction informations on a specific site? ###############
  
-@@ -389,6 +418,7 @@
- syst.attach_lead(lead)
- syst.attach_lead(lead.reversed())
- kwant.plot(syst)
-+save_figure('FAQBA')
- 
- fsyst = syst.finalized()
+ ###### How does Kwant order the propagating modes of a lead? ######
+@@ -348,6 +374,7 @@
+     plt.show()
  
-@@ -401,6 +431,7 @@
- plt.xlabel("momentum [(lattice constant)^-1]")
- plt.ylabel("energy [t]")
- plt.show()
-+save_figure('FAQBB')
+ plot_and_label_modes(flead, E)
++save_figure('FAQPM')
  plt.clf()
  
- wf = kwant.wave_function(fsyst, Ef)
-@@ -474,7 +505,7 @@
+ # More involved example
+@@ -363,6 +390,7 @@
+ flead2 = lead2.finalized()
  
-     # Check that the system looks as intended.
- kwant.plot(syst)
--
-+save_figure('FAQTT')
-     # Finalize the system.
- syst = syst.finalized()
- 
-@@ -488,7 +519,7 @@
- plt.xlabel("momentum [(lattice constant)^-1]")
- plt.ylabel("energy [t]")
- plt.show()
--
-+save_figure('FAQSS')
+ plot_and_label_modes(flead2, 1)
++save_figure('FAQPMC')
  plt.clf()
  
- wf = kwant.wave_function(syst, Ef)
+ 
+@@ -395,7 +423,8 @@
+ 
+ idx = syst.id_by_site[lat(0, 0)]  # look up index of site
+ 
+-print('wavefunction on lat(0, 0): ', wf[idx])
++with open('FAQORD1.txt', 'w') as f:
++    print('wavefunction on lat(0, 0): ', wf[idx], file=f)
+ 
+ lat = kwant.lattice.square(norbs=2)
+ syst = make_system(lat)
+@@ -408,4 +437,5 @@
+ # to degrees of freedom on the same site.
+ wf = wf.reshape(-1, 2)
+ 
+-print('wavefunction on lat(0, 0): ', wf[idx])
++with open('FAQORD2.txt', 'w') as f:
++    print('wavefunction on lat(0, 0): ', wf[idx], file=f)
diff --git a/doc/source/tutorial/FAQ.py b/doc/source/tutorial/FAQ.py
index 0df1b6a9c2e18b413daa69b639b0f312c1387d4d..e8ecbee6d4b0f88834c5999e5ca35d6c4678c5d4 100644
--- a/doc/source/tutorial/FAQ.py
+++ b/doc/source/tutorial/FAQ.py
@@ -1,12 +1,12 @@
+# Frequently Asked Questions
+# ==========================
 
-from cmath import exp
-import numpy as np
 import kwant
-import matplotlib.pyplot
-import matplotlib as mpl
-from matplotlib import pyplot as plt
-matplotlib.rcParams['figure.figsize'] = (3.5,3.5)
+import numpy as np
 import tinyarray
+import matplotlib
+from matplotlib import pyplot as plt
+matplotlib.rcParams['figure.figsize'] = (3.5, 3.5)
 
 
 ######## What is a Site? ##############
@@ -16,8 +16,8 @@ a = 1
 lat = kwant.lattice.square(a)
 syst = kwant.Builder()
 
-syst[lat(1,0)] = 4
-syst[lat(1,1)] = 4
+syst[lat(1, 0)] = 4
+syst[lat(1, 1)] = 4
 
 kwant.plot(syst)
 #HIDDEN_END_FAQ122
@@ -27,47 +27,43 @@ kwant.plot(syst)
 ################# What is a family? a tag? a lattice? ##################
 
 #HIDDEN_BEGIN_FAQ123
-a = 1
-## 2 Monatomic lattices
-lat1 = kwant.lattice.Monatomic( [(1,0) , (0,1)], offset = (0,0)) ## equivalent to kwant.lattice.square()
-lat2 = kwant.lattice.Monatomic( [(1,0) , (0,1)], offset = (0.5,0.5))
+# 2 Monatomic lattices
+primitive_vectors = [(1, 0), (0, 1)]
+lat1 = kwant.lattice.Monatomic(primitive_vectors, offset=(0, 0))  # equivalent to kwant.lattice.square()
+lat2 = kwant.lattice.Monatomic(primitive_vectors, offset=(0.5, 0.5))
 
-## 1 Polyatomic lattice containing two sublattices
-lat3 = kwant.lattice.Polyatomic([(1,0) , (0,1)], [(0,0) , (0.5,0.5)])
+# 1 Polyatomic lattice containing two sublattices
+lat3 = kwant.lattice.Polyatomic([(1, 0) , (0, 1)], [(0, 0) , (0.5, 0.5)])
 subA, subB = lat3.sublattices
 
 
 syst = kwant.Builder()
 
-syst[lat1(0,0)] = 4 ## syst[subA(0,0)] = 4
-syst[lat2(0,0)] = 4 ## syst[subB(0,0)] = 4
+syst[lat1(0, 0)] = 4  # syst[subA(0, 0)] = 4
+syst[lat2(0, 0)] = 4  # syst[subB(0, 0)] = 4
 
 kwant.plot(syst)
 #HIDDEN_END_FAQ123
 
 
-
 ########## What is a hopping? #######################
 
 a = 1
 lat = kwant.lattice.square(a)
 syst = kwant.Builder()
 
-syst[lat(1,0)] = 4
-syst[lat(1,1)] = 4
+syst[lat(1, 0)] = 4
+syst[lat(1, 1)] = 4
 #HIDDEN_BEGIN_FAQ124
-syst[(lat(1,0) , lat(1,1))] = -1
+syst[(lat(1, 0), lat(1, 1))] = -1
 #HIDDEN_END_FAQ124
 
 kwant.plot(syst)
 
-
-
 ########### How to make a hole in a system? ###################"
 
 #HIDDEN_BEGIN_FAQ2
-
-## Definition of the lattice and the system:
+# Define the lattice and the (empty) system
 a = 2
 lat = kwant.lattice.cubic(a)
 syst = kwant.Builder()
@@ -76,47 +72,43 @@ L = 10
 W = 10
 H = 2
 
-## Adding sites to the system:
+# Add sites to the system in a cuboid
 
-syst[ (lat(i,j,k) for i in range(L) for j in range(W) for k in range(H)) ] = 4
+syst[(lat(i, j, k) for i in range(L) for j in range(W) for k in range(H))] = 4
 kwant.plot(syst)
 
+# Delete sites to create a hole
 
-## Deleting sites to create a hole
+def in_hole(site):
+    x, y, z = site.pos / a - (L/2, W/2, H/2)  # position relative to centre
+    return abs(x) < L / 4 and abs(y) < W / 4
+
+for site in filter(in_hole, list(syst.sites())):
+    del syst[site]
 
-for i in range(L):
-    for j in range(W):
-        for k in range(H):
-            x, y, z = lat(i,j,k).pos
-            if ((L-2)*a/4 <= x <= 3*L*a/4) and (W*a/4 <= y <= 3*(W-2)*a/4) and (0 <= z <= (H-1)*a):
-                del syst[lat(i,j,k)]
 kwant.plot(syst)
 #HIDDEN_END_FAQ2
 
 
-
 ################ How can we get access to the sites of our system? ####################
 
 builder = kwant.Builder()
 lat = kwant.lattice.square()
-builder[(lat(i,j) for i in range(3) for j in range(3))] = 4
+builder[(lat(i, j) for i in range(3) for j in range(3))] = 4
 #HIDDEN_BEGIN_FAQ3
-## Before finalizing the system
+# Before finalizing the system
 
-sites = []
-for site in builder.sites():
-    sites.append(site) ## here we choose to add the sites to a list
+sites = list(builder.sites())  # sites() doe *not* return a list
 
 #HIDDEN_END_FAQ3
 #HIDDEN_BEGIN_FAQ7
-## After finalizing the system
-
+# After finalizing the system
 syst = builder.finalized()
-i = syst.id_by_site[lat(0,2)] ## we want the id of the site lat(0,2)
- # syst.sites[i].tag  Returns the tag of lat(0,2)
- # syst.sites[i].pos  Returns the pos of lat(0,2)
+sites = syst.sites  # syst.sites is an actual list
 #HIDDEN_END_FAQ7
-
+#HIDDEN_BEGIN_FAQ72
+i = syst.id_by_site[lat(0, 2)]  # we want the id of the site lat(0, 2)
+#HIDDEN_END_FAQ72
 
 
 ################ How to plot a polyatomic lattice with different colors? ##############"
@@ -125,11 +117,11 @@ i = syst.id_by_site[lat(0,2)] ## we want the id of the site lat(0,2)
 lat = kwant.lattice.kagome()
 syst = kwant.Builder()
 
-a , b, c = lat.sublattices ## The kagome lattice has 3 sublattices
+a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
 #HIDDEN_END_FAQ8
 
 #HIDDEN_BEGIN_FAQ9
-## Plot sites from different families in different colors
+# Plot sites from different families in different colors
 
 def plot_system(syst):
 
@@ -147,14 +139,12 @@ def plot_system(syst):
     kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw)
 
 
-
-
 ## Adding sites and hoppings
 for i in range(4):
     for j in range (4):
-        syst[a(i,j)] = 4 ## red
-        syst[b(i,j)] = 4 ## green
-        syst[c(i,j)] = 4 ## blue
+        syst[a(i, j)] = 4  # red
+        syst[b(i, j)] = 4  # green
+        syst[c(i, j)] = 4  # blue
 
 syst[lat.neighbors()] = -1
 
@@ -164,29 +154,29 @@ plot_system(syst)
 #HIDDEN_END_FAQ9
 
 
-
 ############### How to create every hoppings in a given direction using Hoppingkind? ################
 
-## Monatomic lattice
+# Monatomic lattice
 
 #HIDDEN_BEGIN_FAQ4
 
-## Creation of hopping between neighbors with HoppingKind
+# Create hopping between neighbors with HoppingKind
 a = 1
 syst = kwant.Builder()
 lat = kwant.lattice.square(a)
-syst[ (lat(i,j) for i in range(5) for j in range(5)) ] = 4
+syst[ (lat(i, j) for i in range(5) for j in range(5)) ] = 4
 
-syst[kwant.builder.HoppingKind((1,0), lat)] = -1
+syst[kwant.builder.HoppingKind((1, 0), lat)] = -1
 kwant.plot(syst)
 #HIDDEN_END_FAQ4
 
-## Polyatomic lattice
+# Polyatomic lattice
 
 lat = kwant.lattice.kagome()
 syst = kwant.Builder()
 
-a , b, c = lat.sublattices ## The kagome lattice has 3 sublattices
+a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
+
 
 def plot_system(syst):
 
@@ -198,77 +188,77 @@ def plot_system(syst):
         else:
             return 'green'
 
+    kwant.plot(syst, site_size=0.15, site_lw=0.05, site_color=family_color)
 
-    kwant.plot(syst, site_size= 0.15 , site_lw=0.05, site_color=family_color)
 
 for i in range(4):
     for j in range (4):
-        syst[a(i,j)] = 4 ## red
-        syst[b(i,j)] = 4 ## green
-        syst[c(i,j)] = 4 ## blue
-
+        syst[a(i, j)] = 4  # red
+        syst[b(i, j)] = 4  # green
+        syst[c(i, j)] = 4  # blue
 
 
 #HIDDEN_BEGIN_FAQ13
-syst[kwant.builder.HoppingKind((0,1), b, b)] = -1 # equivalent to syst[kwant.builder.HoppingKind((0,1), b)] = -1
+# equivalent to syst[kwant.builder.HoppingKind((0, 1), b)] = -1
+syst[kwant.builder.HoppingKind((0, 1), b, b)] = -1
 #HIDDEN_END_FAQ13
 plot_system(syst)
-del syst[kwant.builder.HoppingKind((0,1), b, b)] ## delete the hoppings previously created
+# Delete the hoppings previously created
+del syst[kwant.builder.HoppingKind((0, 1), b, b)]
 #HIDDEN_BEGIN_FAQ14
-syst[kwant.builder.HoppingKind((0,0), a, b)] = -1
-syst[kwant.builder.HoppingKind((0,0), a, c)] = -1
-syst[kwant.builder.HoppingKind((0,0), c, b)] = -1
+syst[kwant.builder.HoppingKind((0, 0), a, b)] = -1
+syst[kwant.builder.HoppingKind((0, 0), a, c)] = -1
+syst[kwant.builder.HoppingKind((0, 0), c, b)] = -1
 #HIDDEN_END_FAQ14
 plot_system(syst)
 
 
-
 ########## How to create the hoppings between adjacent sites? ################
 
-## Monatomic lattice
+# Monatomic lattice
 
 #HIDDEN_BEGIN_FAQ5
 
-## Creation of hoppings with lat.neighbors()
+# Create hoppings with lat.neighbors()
 syst = kwant.Builder()
 lat = kwant.lattice.square()
-syst[ (lat(i,j) for i in range(3) for j in range(3)) ] = 4
+syst[(lat(i, j) for i in range(3) for j in range(3))] = 4
 
-syst[lat.neighbors()] = -1 ## equivalent to lat.neighbors(1)
+syst[lat.neighbors()] = -1  # Equivalent to lat.neighbors(1)
 kwant.plot(syst)
 
-del syst[lat.neighbors()] ## deletes every hoppings previously created to add new one
+del syst[lat.neighbors()]  # Delete all nearest-neighbor hoppings
 syst[lat.neighbors(2)] = -1
 
 kwant.plot(syst)
 #HIDDEN_END_FAQ5
 
-## Polyatomic lattice
+# Polyatomic lattice
 
 #HIDDEN_BEGIN_FAQ6
 
-## Hoppings using .neighbors()
+# Hoppings using .neighbors()
 #HIDDEN_BEGIN_FAQ10
-## Creation of the system
+# Create the system
 lat = kwant.lattice.kagome()
 syst = kwant.Builder()
-a , b, c = lat.sublattices ## The kagome lattice has 3 sublattices
+a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
+
 for i in range(4):
     for j in range (4):
-        syst[a(i,j)] = 4 ## red
-        syst[b(i,j)] = 4 ## green
-        syst[c(i,j)] = 4 ## blue
+        syst[a(i, j)] = 4  # red
+        syst[b(i, j)] = 4  # green
+        syst[c(i, j)] = 4  # blue
 
-## .neighbors()
 syst[lat.neighbors()] = -1
 #HIDDEN_END_FAQ10
 plot_system(syst)
-del syst[lat.neighbors()] ## delete the hoppings previously created
+del syst[lat.neighbors()]  # Delete the hoppings previously created
 #HIDDEN_BEGIN_FAQ11
 syst[a.neighbors()] = -1
 #HIDDEN_END_FAQ11
 plot_system(syst)
-del syst[a.neighbors()] ## deletes every hoppings previously created to add new one
+del syst[a.neighbors()]  # Delete the hoppings previously created
 
 
 #HIDDEN_BEGIN_FAQ12A
@@ -278,12 +268,9 @@ plot_system(syst)
 del syst[lat.neighbors(2)]
 
 
-
 ##### How to create a lead with a lattice different from the scattering region? ##########
 
-
-#HIDDEN_BEGIN_FAQAA
-## Plot sites from different families in different colors
+# Plot sites from different families in different colors
 
 def plot_system(syst):
 
@@ -295,94 +282,54 @@ def plot_system(syst):
         else:
             return 'green'
 
-
     kwant.plot(syst, site_lw=0.1, site_color=family_color)
 
-#HIDDEN_END_FAQAA
-#HIDDEN_BEGIN_FAQAB
-## Defining the scattering Region
-a = 2
-lat = kwant.lattice.honeycomb(a)
-syst = kwant.Builder()
 
+#HIDDEN_BEGIN_FAQAA
+# Define the scattering Region
 L = 5
 W = 5
 
+lat = kwant.lattice.honeycomb()
 subA, subB = lat.sublattices
 
-## Adding sites to the system:
-
-syst[ (subA(i,j) for i in range(L) for j in range(W)) ] = 4
-syst[ (subB(i,j) for i in range(L) for j in range(W)) ] = 4
+syst = kwant.Builder()
+syst[(subA(i, j) for i in range(L) for j in range(W))] = 4
+syst[(subB(i, j) for i in range(L) for j in range(W))] = 4
 syst[lat.neighbors()] = -1
+#HIDDEN_END_FAQAA
 plot_system(syst)
-#HIDDEN_END_FAQAB
-#HIDDEN_BEGIN_FAQAC
-## We manually add sites of the same lead lattice
 
+#HIDDEN_BEGIN_FAQAB
+# Create a lead
+lat_lead = kwant.lattice.square()
+sym_lead1 = kwant.TranslationalSymmetry((0, 1))
 
-lat2 = kwant.lattice.square(a)
-
-def shapetop(pos):
-    x, y = pos
-    return ( 4 <= x <= 12 ) and ( 8 < y <= 10 )
-
-def shapebot(pos):
-    x, y = pos
-    return ( 0 <= x <= 8 ) and ( -2 <= y < 0 )
-
-
-syst[lat2.shape(shapetop, (4,10))] = 4
-syst[lat2.shape(shapebot, (0,-2))] = 4
-
-syst[lat2.neighbors()] = -1
-syst[((lat2(i,-1), subA(i,0)) for i in range(5))] = -1
-syst[((lat2(i+2,5), subB(i,4)) for i in range(5))] = -1
-plot_system(syst)
-#HIDDEN_END_FAQAC
-#HIDDEN_BEGIN_FAQAD
-## Creation of the top lead
-
-lat_lead = kwant.lattice.square(a)
-sym_lead1 = kwant.TranslationalSymmetry((0,a))
 lead1 = kwant.Builder(sym_lead1)
+lead1[(lat_lead(i, 0) for i in range(2, 7))] = 4
+lead1[lat_lead.neighbors()] = -1
+#HIDDEN_END_FAQAB
+plot_system(lead1)
 
+#HIDDEN_BEGIN_FAQAC
+syst[(lat_lead(i, 5) for i in range(2, 7))] = 4
+syst[lat_lead.neighbors()] = -1
 
-def lead_shape_top(pos): ## Shape of the lead
-    (x, y) = pos
-    return (4 <= x <= 12)
-
-
-lead1[lat_lead.shape(lead_shape_top, (4,12))] = 4
-lead1[lat_lead.neighbors()] = -1
+# Manually attach sites from graphene to square lattice
+syst[((lat_lead(i+2, 5), subB(i, 4)) for i in range(5))] = -1
+#HIDDEN_END_FAQAC
+plot_system(syst)
 
+#HIDDEN_BEGIN_FAQAD
 syst.attach_lead(lead1)
-plot_system(syst)
 #HIDDEN_END_FAQAD
-#HIDDEN_BEGIN_FAQAE
-## Creation of the bottom lead
-
-sym_lead2 = kwant.TranslationalSymmetry((0,-a))
-lead2 = kwant.Builder(sym_lead2)
-
-def lead_shape_bot(pos): ## Shape of the lead
-    (x, y) = pos
-    return (0 <= x <= 8)
-
-lead2[lat_lead.shape(lead_shape_bot, (0,-4))] = 4
-lead2[lat_lead.neighbors()] = -1
-
-syst.attach_lead(lead2)
-
 plot_system(syst)
-#HIDDEN_END_FAQAE
-
 
 
 ############# How to cut a finite system out of a system with translationnal symmetries? ###########
 
 #HIDDEN_BEGIN_FAQccc
-# Create  3d model.
+# Create 3d model.
 cubic = kwant.lattice.cubic()
 sym_3d = kwant.TranslationalSymmetry([1, 0, 0], [0, 1, 0], [0, 0, 1])
 model = kwant.Builder(sym_3d)
@@ -398,9 +345,8 @@ def cuboid_shape(site):
 
 cuboid = kwant.Builder()
 cuboid.fill(model, cuboid_shape, (0, 0, 0));
-
-kwant.plot(cuboid);
 #HIDDEN_END_FAQddd
+kwant.plot(cuboid);
 
 #HIDDEN_BEGIN_FAQeee
 # Build electrode (black).
@@ -409,164 +355,110 @@ def electrode_shape(site):
     return y**2 + z**2 < 2.3**2
 
 electrode = kwant.Builder(kwant.TranslationalSymmetry([1, 0, 0]))
-electrode.fill(model, electrode_shape, (0, 5, 2)); ## lead
+electrode.fill(model, electrode_shape, (0, 5, 2))  # lead
 
-cuboid.fill(electrode, lambda s: abs(s.pos[0]) < 7, (0, 5, 4)); ## scattering region
+# Scattering region
+cuboid.fill(electrode, lambda s: abs(s.pos[0]) < 7, (0, 5, 4))
 
 cuboid.attach_lead(electrode)
-kwant.plot(cuboid);
 #HIDDEN_END_FAQeee
+kwant.plot(cuboid);
 
 
-###### How to extract the wavefunction informations on a specific site? ###############
+###### How does Kwant order the propagating modes of a lead? ######
 
-#HIDDEN_BEGIN_FAQAF
-## Creation of the system
-a = 2
-lat = kwant.lattice.square(a)
-syst = kwant.Builder()
-
-syst[((lat(i,j) for i in range(5) for j in range(3)))] = 4
-syst[lat.neighbors()] = -1
-kwant.plot(syst)
+#HIDDEN_BEGIN_PM
+lat = kwant.lattice.square()
 
-sym_lead = kwant.TranslationalSymmetry((-a,0))
-lead = kwant.Builder(sym_lead)
-lead[(lat(0,i) for i in range(3))] = 4
+lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+lead[(lat(0, i) for i in range(3))] = 4
 lead[lat.neighbors()] = -1
-syst.attach_lead(lead)
-syst.attach_lead(lead.reversed())
-kwant.plot(syst)
 
-fsyst = syst.finalized()
-
-## Plot the different modes
-Ef = 3.0
-lead = lead.finalized()
-kwant.plotter.bands(lead, show=False)
-kwant.plotter.bands
-plt.plot([-3,3],[Ef,Ef], 'r--')
-plt.xlabel("momentum [(lattice constant)^-1]")
-plt.ylabel("energy [t]")
-plt.show()
-#HIDDEN_END_FAQAF
+flead = lead.finalized()
+
+E = 2.5
+prop_modes, _ = flead.modes(energy=E)
+#HIDDEN_END_PM
+
+def plot_and_label_modes(lead, E):
+    # Plot the different modes
+    pmodes, _ = lead.modes(energy=E)
+    kwant.plotter.bands(lead, show=False)
+    for i, k in enumerate(pmodes.momenta):
+        plt.plot(k, E, 'ko')
+        plt.annotate(str(i), xy=(k, E), xytext=(-5, 8),
+                     textcoords='offset points',
+                     bbox=dict(boxstyle='round,pad=0.1',fc='white', alpha=0.7))
+    plt.plot([-3, 3], [E, E], 'r--')
+    plt.ylim(E-1, E+1)
+    plt.xlim(-2, 2)
+    plt.xlabel("momentum")
+    plt.ylabel("energy")
+    plt.show()
+
+plot_and_label_modes(flead, E)
 plt.clf()
 
-#HIDDEN_BEGIN_FAQAG
-wf = kwant.wave_function(fsyst, Ef)
-wf_left_lead = wf(0) ## Wave function for the first lead (left)
-#HIDDEN_END_FAQAG
-
+# More involved example
 
-#HIDDEN_BEGIN_FAQAH
+s0 = np.eye(2)
+sz = np.array([[1, 0], [0, -1]])
 
-wf_left_0 = wf_left_lead[0] ## We choose the the mode with the highest k (mode in blue)
+lead2 = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
 
+lead2[(lat(0, i) for i in range(2))] = np.diag([1.8, -1])
+lead2[lat.neighbors()] = -1 * sz
 
-tag = lat.closest((6,2)) ## returns the tag  of the site from lat based on the position
-
-i = fsyst.id_by_site[lat(*tag)] ## Returns the index in the low_level system based on the tag
-
-wf_site = wf_left_0[i] ## Returns the wave function on the site
-
-#HIDDEN_END_FAQAH
+flead2 = lead2.finalized()
 
+plot_and_label_modes(flead2, 1)
+plt.clf()
 
-#HIDDEN_BEGIN_FAQAO
-tau_x = tinyarray.array([[0, 1], [1, 0]])
-tau_y = tinyarray.array([[0, -1j], [1j, 0]])
-tau_z = tinyarray.array([[1, 0], [0, -1]])
 
-lat = kwant.lattice.square(norbs=2)
+###### How does Kwant order components of an individual wavefunction? ######
 
+def circle(R):
+    return lambda r: np.linalg.norm(r) < R
 
-	### Creation of the Builder ###
 
-def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
-                mu=0.4, Delta=0.1, Deltapos=4, t=1.0, phs=True):
-    # Start with an empty tight-binding system. On each site, there
-    # are now electron and hole orbitals, so we must specify the
-    # number of orbitals per site. The orbital structure is the same
-    # as in the Hamiltonian.
+def make_system(lat):
+    norbs = lat.norbs
     syst = kwant.Builder()
+    syst[lat.shape(circle(3), (0, 0))] = 4 * np.eye(norbs)
+    syst[lat.neighbors()] = -1 * np.eye(norbs)
 
-    #### Define the scattering region. ####
-    # The superconducting order parameter couples electron and hole orbitals
-    # on each site, and hence enters as an onsite potential.
-    # The pairing is only included beyond the point 'Deltapos' in the scattering region.
-    syst[(lat(x, y) for x in range(Deltapos) for y in range(W))] = (4 * t - mu) * tau_z
-    syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
-
-    # The tunnel barrier
-    syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1])
-         for y in range(W))] = (4 * t + barrier - mu) * tau_z
-
-    # Hoppings
-    syst[lat.neighbors()] = -t * tau_z
-    #### Define the leads. ####
-    # Left lead - normal, so the order parameter is zero.
-    sym_left = kwant.TranslationalSymmetry((-a, 0))
-    # Specify the conservation law used to treat electrons and holes separately.
-    # We only do this in the left lead, where the pairing is zero.
-    lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y)
-    lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z
-    lead0[lat.neighbors()] = -t * tau_z
-    # Right lead - superconducting, so the order parameter is included.
-    sym_right = kwant.TranslationalSymmetry((a, 0))
-    lead1 = kwant.Builder(sym_right)
-    lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
-    lead1[lat.neighbors()] = -t * tau_z
-
-    #### Attach the leads and return the system. ####
-    syst.attach_lead(lead0)
-    syst.attach_lead(lead1)
-
-    return syst
-
-
-syst = make_system(W=10)
-
-    # Check that the system looks as intended.
-kwant.plot(syst)
-
-    # Finalize the system.
-syst = syst.finalized()
-
-	#Plot the mods of the lead
-lead = syst.leads[0]
-
-Ef=0.5
-kwant.plotter.bands(lead, show=False)
-kwant.plotter.bands
-plt.plot([-3,3],[Ef,Ef], 'r--')
-plt.xlabel("momentum [(lattice constant)^-1]")
-plt.ylabel("energy [t]")
-plt.show()
-
-plt.clf()
+    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+    lead[(lat(0, i) for i in range(-1, 2))] = 4 * np.eye(norbs)
+    lead[lat.neighbors()] = -1 * np.eye(norbs)
 
-wf = kwant.wave_function(syst, Ef)
-wf_left_lead = wf(0) # Wave function for the first lead (left)
+    syst.attach_lead(lead)
+    syst.attach_lead(lead.reversed())
 
+    return syst.finalized()
 
 
-#HIDDEN_END_FAQAO
+#HIDDEN_BEGIN_ORD1
+lat = kwant.lattice.square(norbs=1)
+syst = make_system(lat)
+scattering_states = kwant.wave_function(syst, energy=1)
+wf = scattering_states(0)[0]  # scattering state from lead 0 incoming in mode 0
 
+idx = syst.id_by_site[lat(0, 0)]  # look up index of site
 
-#HIDDEN_BEGIN_FAQAP
-nb_degrees = 2 ## 2 degrees of freedom
+print('wavefunction on lat(0, 0): ', wf[idx])
+#HIDDEN_END_ORD1
 
-wf_left_0 = wf_left_lead[0] ## We choose the the mode with the highest k (mode in blue)
-
-
-tag = lat.closest((6,2)) ## returns the tag  of the site from lat based on the position
-
-i = syst.id_by_site[lat(*tag)] ## Returns the index in the low_level system based on the tag
-
-tab_wf = []
+#HIDDEN_BEGIN_ORD2
+lat = kwant.lattice.square(norbs=2)
+syst = make_system(lat)
+scattering_states = kwant.wave_function(syst, energy=1)
+wf = scattering_states(0)[0]  # scattering state from lead 0 incoming in mode 0
 
-for k in range(nb_degrees) : ## loop over the number of degrees of freedom
-	tab_wf.append(wf_left_0[nb_degrees * i +k]) # different states of the given site
+idx = syst.id_by_site[lat(0, 0)]  # look up index of site
 
+# Group consecutive degrees of freedom from 'wf' together; these correspond
+# to degrees of freedom on the same site.
+wf = wf.reshape(-1, 2)
 
-#HIDDEN_END_FAQAP
+print('wavefunction on lat(0, 0): ', wf[idx])
+#HIDDEN_END_ORD2
diff --git a/doc/source/tutorial/FAQ.rst b/doc/source/tutorial/FAQ.rst
index 48ec9abadb011ca18846b6c0e0332a178e3f9614..58dbc096fa6c03f8531cc60e3d59fd8138d8c248 100644
--- a/doc/source/tutorial/FAQ.rst
+++ b/doc/source/tutorial/FAQ.rst
@@ -1,17 +1,37 @@
----------------------------------------------------------------
-Frequently Asked Questions :
----------------------------------------------------------------
+--------------------------
+Frequently Asked Questions
+--------------------------
+It is important to read the tutorials before looking at the questions. This FAQ
+is aimed to add complementary explainations that are not in the tutorials. The `Kwant paper <https://downloads.kwant-project.org/doc/kwant-paper.pdf>`_ also digs deeper into Kwant's structure.
 
-It is important to read the tutorials before looking at the questions. This FAQ is made to add complementary explanations that are not in the tutorials.
 
+What is a system, and what is a builder?
+========================================
+A Kwant system represents a tight-binding model. It contains a graph of edges
+and vertices that are assigned values, and which is used to construct
+the Hamiltonian for the model being simulated.
 
-What is a site?
-=================
+In Kwant the specification of the tight-binding model is separated from the use
+of the model in numerical calculations. The `~kwant.builder.Builder` is used
+when specifying/constructing the model, then the
+`~kwant.builder.Builder.finalize` method is called, which produces a so-called
+low-level `~kwant.system.System` that can be used by Kwant's solvers.
+
+In the documentation and in mailing list discussions, the term general term
+"system" can be used to refer either to a ``Builder`` or to a low-level
+``System``, and the context will determine which specific class is being
+referred to. The terms "builder" and "low-level system" (or "finalized system")
+refer respectively to ``Builder`` and ``System``.
 
-Kwant is based on the tight-binding model which is composed of edges and vertices.
-Site objects are Kwant’s abstraction for the vertices. The sites have two attributes: the **family** and the **tag** .
 
-For example :
+What is a site?
+===============
+Kwant is a tool for working with tight-binding models, which can be viewed as a
+graph composed of edges and vertices.  Site objects are Kwant’s abstraction for
+the vertices.  Sites have two attributes: a **family** and a **tag** .  The
+combination of family and tag uniquely define a site.
+
+For example let us create an empty tight binding system and add two sites:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQ122
@@ -19,18 +39,45 @@ For example :
 
 .. image:: ../images/FAQ122.*
 
-As we can see in the figure, we added 2 sites to the system. Each site has a different tag ( ``(1,0)`` and ``(1,1)`` ) but they have the same family ``lat``
+In the above snippet we added 2 sites: ``lat(1 ,0)`` and ``lat(0 , 1)``. Both
+of these sites belong to the same family, ``lat``, but have different tags:
+``(1, 0)`` and ``(0, 1)`` respectively.
 
-What is a lattice?
-====================================
 
-The lattice contains the spatial repartition of the sites in the system. There are two kind of lattices in Kwant :
-	- The monatomic lattices
-	- The polyatomic lattices
+What is a site family, and what is a tag?
+=========================================
+a site family "groups" related sites together, and a tag serves as a unique
+identifier for a site within a given family.
+
+In the previous example we saw a family that was suggestively called ``lat``,
+which had sites whose tags were pairs of integers. In this specific example
+the site family also happens to be a regular Bravais lattice, and the tags take
+on the meaning of lattice coordinates for a site on this lattice.
 
+The concept of families and tags is, however, more general. For example, one
+could implement a mesh that can be locally refined in certain areas, by having
+a family where sites belong to a `quadtree
+<https://en.wikipedia.org/wiki/Quadtree>`_, or an amorphous blob where sites
+are tagged by letters of the alphabet.
 
-The monatomic class represents a bravais lattice with a single site in the basis, it contains an origin and multiple vectors which define a frame of reference. The polyatomic lattice has N sites per basis which correspond to N monatomic sublattices.
 
+What is a lattice?
+==================
+Kwant allows users to define and use Bravais lattices for dealing with
+collections of regularly placed sites. They know about things like which sites
+on the lattice are neighbors, and how to fill a region of realspace with sites.
+``Monatomic`` lattices have a single site in their basis, while ``Polyatomic``
+have more than one site in their basis.
+
+``Monatomic`` lattices in Kwant *are also site families*, with sites that are
+tagged with tuples of integers: the site's coordinates in the basis of
+primitive vectors of the lattice. ``Polyatomic`` lattices, however, are *not*
+site families, as lattice coordinates are not enough information to uniquely
+identify a site if there is more than one site in the basis. ``Polyatomic``
+lattices do, however, have an attribute ``sublattices`` that is a list of
+``Monatomic`` lattices that together make up the whole ``Polyatomic`` lattice.
+
+For example:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQ123
@@ -38,91 +85,70 @@ The monatomic class represents a bravais lattice with a single site in the basis
 
 .. image:: ../images/FAQ123.*
 
-As we can see above, we created 2 monatomic lattices,  ``(1,0)`` and ``(0,1)`` are the primitive vectors and ``(0,0)`` and ``(0.5,0.5)`` are the origins of the two lattices. We can create a polyatomic lattice with the same primitive vectors and 2 sites in the basis. It leads to the same result. The two sublattices are equivalent to the two monatomic lattices.
+Above, we created 2 ``Monatomic`` lattices (``lat1`` and ``lat2``).  ``(1, 0)``
+and ``(0, 1)`` are the primitive vectors and ``(0, 0)`` and ``(0.5, 0.5)`` are
+the origins of the two lattices. Next we create a ``Polyatomic`` lattice with the
+same primitive vectors and 2 sites in the basis.
 
-What is a family ? a tag?
-============================
+The two sublattices are equivalent to the two monatomic lattices that we
+created previously. Because a ``Polyatomic`` lattice knows about its
+sublattices, we can do things like calculate neighboring sites, even between
+sublattices, which would not be possible with the two separate ``Monatomic``
+lattices.
 
-The concept of lattices and families are very interrelated. The family represents a sublattice and defines the tag.
-In the monatomic case, the family has the same name as the the lattice.
+the `kwant.lattice` module also defines several functions, such as
+`~kwant.lattice.square` and `~kwant.lattice.honeycomb`, which serve as a
+convenience for creating lattices of common types, without having to
+explicitly specify all of the lattice vectors and basis vectors.
 
-In the polyatomic case, there are multiple families which represent the different monatomic sublattices.
 
-If we take the example from above, there are 4 families: ``subA , subB , lat1, lat2``. However, ``lat3`` doesn't have a family because it is polyatomic.
-
-The tag represents the coordinates of a site in his family. In the common case where the site family is a regular lattice, the site tag is simply given by its integer lattice coordinates.
-
-There are multiple :doc:`predefined lattices <../reference/kwant.lattice>` implemented in Kwant that we can use, it is simple to see how much sublattices there are by using ``lat.sublattices`` and look at what it returns.
-
-
-What is a hopping?
-====================
-
-If we take the example of `What is a site?`_ we can simply add a hopping with :
+How do I plot a polyatomic lattice with different colors?
+=========================================================
+In the following example we shall use a kagome lattice, which has 3 sublattices.
 
 .. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQ124
-    :end-before: #HIDDEN_END_FAQ124
-
-.. image:: ../images/FAQ124.*
-
-As we can see, a hopping is a tupple of site, it represents an edge of the tight-binding system.
-
-How to make a hole in a system?
-=================================
-To make a hole in the system, we use ``del syst[key]`` , we can either use the real-space coordinates with ``.pos`` or use the lattice coordinates ``tag`` . Here, we use the real-space position to delete sites in the middle.
-
-.. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQ2
-    :end-before: #HIDDEN_END_FAQ2
-
-.. image:: ../images/FAQ1.*
-.. image:: ../images/FAQ2.*
-
-
-Note that we can make the hole after creating the hoppings, then the hoppings to the sites that are deleted are removed as well.
-
-How can we get access to the sites of our system?
-===================================================
-
-In kwant, the tight-binding model is split in two steps. First, we have a Builder where it is easy to add, remove or modify sites and hoppings. Then, we have a System which is optimized for computation. The separation is made with the finalization. The way we can have access to the list of sites is different before and after the finalization.
+    :start-after: #HIDDEN_BEGIN_FAQ8
+    :end-before: #HIDDEN_END_FAQ8
 
+As we can see below, we create a new plotting function that assigns a color for each family, and a different size for the hoppings depending on the family of the two sites. Finally we add sites and hoppings to our system and plot it with the new function.
 
 .. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQ3
-    :end-before: #HIDDEN_END_FAQ3
+    :start-after: #HIDDEN_BEGIN_FAQ9
+    :end-before: #HIDDEN_END_FAQ9
 
-Before finalizing, we get access to the sites with : ``syst.sites()`` .
+.. image:: ../images/FAQ6B.*
 
-.. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQ7
-    :end-before: #HIDDEN_END_FAQ7
 
-After finalizing the system, the order of the sites changes, so we need ``id_by_site`` to find the site's id in the finalized system. We can also use ``lat.closest(pos)`` to get the tag of a site at a given position.
+What is a hopping?
+==================
+A hopping is simply a pair of sites, which defines an edge of the graph
+that makes up our tight-binding model.
 
-How to plot a polyatomic lattice with different colors?
-===============================================================================
-We take a kagome lattice as example here, it has 3 sublattices.
+Starting from the example code from `What is a site?`_, we can add a hopping
+to our system in the following way:
 
 .. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQ8
-    :end-before: #HIDDEN_END_FAQ8
-
-As we can see below, we create a new plotting function that assigns a color for each family, and a different size for the hoppings depending on the family of the two sites. Finally we add sites and hoppings to our system and plot it with the new function.
+    :start-after: #HIDDEN_BEGIN_FAQ124
+    :end-before: #HIDDEN_END_FAQ124
 
-.. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQ9
-    :end-before: #HIDDEN_END_FAQ9
+.. image:: ../images/FAQ124.*
 
-.. image:: ../images/FAQ6B.*
+Visually, we represent a hopping as a line that joins two sites.
 
-How to create every hoppings in a given direction using ``Hoppingkind``?
-==========================================================================
+Whenever we add a hopping ``(site_a, site_b)`` to a system and assign it a
+value ``v``, implicitly the hopping ``(site_b, site_a)`` is also added, with
+value the Hermitian conjugate of ``v``.
 
 
-If the lattice is monatomic:
+How do I create all hoppings in a given direction?
+==================================================
+This can be obtained using a `~kwant.builder.HoppingKind`. A ``HoppingKind`` is
+a way of specifying all hoppings of a particular "kind", between two site
+families. For example ``HoppingKind((1, 0), lat_a, lat_b)`` represents all
+hoppings of the form ``(lat_a(x + (1, 0)), lat_b(x))``, where ``x`` is a tag
+(here, a pair of integers).
 
-``HoppingKind`` can be used to easily create a hopping in a specific direction.
+The following example shows how this can be used:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQ4
@@ -130,11 +156,9 @@ If the lattice is monatomic:
 
 .. image:: ../images/FAQ3.*
 
-We create hoppings between the sites of the lattice ``lat`` . ``(1,0)`` represents the direction of the hoppings we want to assign.
-
-If the lattice is polyatomic:
-
-``Hoppingkind`` can only be used in the monatomic sublattices. As we can see below, it can create hoppings between different sublattices with a given direction.
+Note that ``HoppingKind`` only works with site families so you cannot use
+them directly with ``Polyatomic`` lattices; you have to explicitly specify
+the sublattices when creating a ``HoppingKind``:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQ13
@@ -153,15 +177,11 @@ Here, we create hoppings between the sites of the same lattice coordinates but f
 .. image:: ../images/FAQ11.*
 
 
-
-How to create the hoppings between adjacent sites?
+How do I create the hoppings between adjacent sites?
 ====================================================
-
-``lat.neighbors(n)`` returns ``HoppingKind`` in the directions of the n-nearest neighbors.
-
-In the monatomic case:
-
-``lat.neighbors(n)`` creates the hoppings between the n-nearest neighbors.
+``Polyatomic`` and ``Monatomic`` lattices have a method `~kwant.lattice.Polyatomic.neighbors`
+that returns a (or several) ``HoppingKind`` that connects sites with their
+(n-nearest) neighors:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQ5
@@ -172,47 +192,96 @@ In the monatomic case:
 
 As we can see in the figure above, ``lat.neighbors()`` (on the left) returns the hoppings between the first nearest neighbors and ``lat.neighbors(2)`` (on the right) returns the hoppings between the second nearest neighbors.
 
-In the polyatomic case:
-
-It is possible to use ``.neighbors()`` with the lattice or the sublattices as explained below.
+When using a ``Polyatomic`` lattice ``neighbors()`` knows about the different
+sublattices:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQ10
     :end-before: #HIDDEN_END_FAQ10
 
-It creates hoppings between the first nearest neighbors of every sites of the system.
 
 .. image:: ../images/FAQ7.*
 
+However, if we use the ``neighbors()`` method of a single sublattice, we will
+only get the neighbors *on that sublattice*:
+
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQ11
     :end-before: #HIDDEN_END_FAQ11
 
-It creates hoppings between the first nearest neighbors of the sites of the sublattice a (in red).
-
 .. image:: ../images/FAQ8.*
 
+Note in the above that there are *only* hoppings between the blue sites. This
+is an artifact of the visualisation: the red and green sites just happen to lie
+in the path of the hoppings, but are not connected by them.
+
+
+How do I make a hole in a system?
+=================================
+To make a hole in the system, we use ``del syst[site]``. In the following
+example we remove all sites inside some "hole" region:
+
 .. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQ12A
-    :end-before: #HIDDEN_END_FAQ12A
+    :start-after: #HIDDEN_BEGIN_FAQ2
+    :end-before: #HIDDEN_END_FAQ2
 
-It creates hoppings between the second nearest neighbors of every sites of the system.
+.. image:: ../images/FAQ1.*
+.. image:: ../images/FAQ2.*
 
-.. image:: ../images/FAQ9.*
+``del syst[site]`` also works after hoppings have been added to the system.
+If a site is deleted, then all the hoppings to/from that site are also deleted.
 
 
+How can I get access to a system's sites?
+=========================================
+The ways of accessing system sites is slightly different depending on whether
+we are talking about a ``Builder`` or ``System`` (see `What is a system, and
+what is a builder?`_ if you do not know the difference).
 
-How to create a lead with a lattice different from the scattering region?
+We can access the sites of a ``Builder`` by using its `~kwant.builder.Builder.sites` method:
+
+.. literalinclude:: FAQ.py
+    :start-after: #HIDDEN_BEGIN_FAQ3
+    :end-before: #HIDDEN_END_FAQ3
+
+The ``sites()`` method returns an *iterator* over the system sites, and in the
+above example we create a list from the contents of this iterator, which
+contains all the sites. At this stage the ordering of sites is not fixed, so if
+you add more sites to the ``Builder`` and call ``sites()`` again, the sites may
+well be returned in a different order.
+
+After finalization, when we are dealing with a ``System``, the sites themselves
+are stored in a list, which can be accessed via the ``sites`` attribute:
+
+.. literalinclude:: FAQ.py
+    :start-after: #HIDDEN_BEGIN_FAQ7
+    :end-before: #HIDDEN_END_FAQ7
+
+The order of sites in a ``System`` is fixed, and also defines the ordering of
+the system Hamiltonian, system wavefunctions etc. (see `How does Kwant order components of an individual wavefunction?`_ for details).
+
+``System`` also contains the inverse mapping, ``id_by_site`` which gives us
+the index of a given site within the system:
+
+.. literalinclude:: FAQ.py
+    :start-after: #HIDDEN_BEGIN_FAQ72
+    :end-before: #HIDDEN_END_FAQ72
+
+
+How do I create a lead with a lattice different from the scattering region?
 ===========================================================================
+Let us take the example of a system containing sites from a honeycomb lattice,
+which we want to connect to leads that contain sites from a square lattice.
 
-First of all, we need to plot the sites in different colors depending on the sites families.
+First we construct the central system:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQAA
     :end-before: #HIDDEN_END_FAQAA
 
+.. image:: ../images/FAQAA.*
 
-Then, we create the scattering region, here we take the example of the honeycomb lattice.
+and the lead:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQAB
@@ -220,7 +289,12 @@ Then, we create the scattering region, here we take the example of the honeycomb
 
 .. image:: ../images/FAQAB.*
 
-We now have the scattering region, however, we can't just make the lead and attach it if it doesn't have the same lattice. We want our lead lattice to be square, so we need to manually attach some sites from the same lead lattice to the scattering region.
+We cannot simply use `~kwant.builder.Builder.attach_lead` to attach this lead to the
+system with the honeycomb lattice because Kwant does not know how sites from
+these two lattices should be connected.
+
+We must first add a layer of sites from the square lattice to the system and manually
+add the hoppings from these sites to the sites from the honeycomb lattice:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQAC
@@ -228,7 +302,7 @@ We now have the scattering region, however, we can't just make the lead and atta
 
 .. image:: ../images/FAQAC.*
 
-Now that we have those sites, we can create the leads and attach it. We first begin with the top lead, considering the dimensions of the scattering region on the figure above, we create a shape for the lead and attach it.
+``attach_lead()`` will now be able to attach the lead:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQAD
@@ -236,33 +310,25 @@ Now that we have those sites, we can create the leads and attach it. We first be
 
 .. image:: ../images/FAQAD.*
 
-We now do the same for the bottom lead.
-
-.. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQAE
-    :end-before: #HIDDEN_END_FAQAE
-
-.. image:: ../images/FAQAE.*
-
-
 
-How to cut a finite system out of a system with translationnal symmetries?
-============================================================================
-
-This can be achieved using the :doc:`fill()  <../reference/generated/kwant.builder.Builder>` method to fill a Builder with a Builder with higher symmetry.
-
-An example of fill is given :doc:`here  <../tutorial/discretize>` for an undisctretised hamiltonian, but we want to see how it works with a discretised system.
-
-When using the fill method, we need two systems, the template and the target (final sytem). The template is a part of the system which is periodically repeated in the desired shape to create the final system, so it needs to have translationnal symmetries.
+How do I cut a finite system out of a system with translational symmetries?
+===========================================================================
+This can be achieved using the `~kwant.builder.Builder.fill` method to fill a
+``Builder`` with a ``Builder`` with higher symmetry.
 
+When using the ``fill()`` method, we need two systems: the template and the
+target. The template is a ``Builder`` with some translational symmetry that
+will be repeated in the desired shape to create the final system.
 
-We first create the system template with its symmetries.
+For example, say we want to create a simple model on a cubic lattice:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQccc
     :end-before: #HIDDEN_END_FAQccc
 
-We now define the shape of the scattering region and use the fill method to create the first part of the scattering region
+We have now created our "template" ``Builder`` which has 3 translational
+symmetries. Next we will fill a system with no translational symmetries with
+sites and hoppings from the template inside a cuboid:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQddd
@@ -270,7 +336,9 @@ We now define the shape of the scattering region and use the fill method to crea
 
 .. image:: ../images/FAQaaa.*
 
-Finally, we create the template that will used to create the second part of the scattering region and the lead.
+We can then use the original template to create a lead, which has 1 translational
+symmetry. We can then use this lead as a template to fill another section of
+the system with a cylinder of sites and hoppings:
 
 .. literalinclude:: FAQ.py
     :start-after: #HIDDEN_BEGIN_FAQeee
@@ -278,58 +346,84 @@ Finally, we create the template that will used to create the second part of the
 
 .. image:: ../images/FAQbbb.*
 
-How are the physical quantities ordered within the output of the functions?
-===========================================================================
-
-We take the example of the ``wave_function()`` here. We have only one degree of freedom here, see `How are degrees of freedom ordered ?`_  for further explanations in case of multiple one.
 
-So first, we define a system with a square lattice and 2 leads. We also plot the Fermi energy and the modes to see how much modes we can get with a given energy.
+How does Kwant order the propagating modes of a lead?
+=====================================================
+A very useful feature of kwant is to calculate the transverse wavefunctions of
+propagating modes in a system with 1 translational symmetry.  This can be
+achieved with the `~kwant.system.InfiniteSystem.modes` method, which returns a
+pair of objects, the first of which contains the propagating modes of the
+system in a `~kwant.physics.PropagatingModes` object:
 
 .. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQAF
-    :end-before: #HIDDEN_END_FAQAF
+    :start-after: #HIDDEN_BEGIN_PM
+    :end-before: #HIDDEN_END_PM
 
-.. image:: ../images/FAQBA.*
-.. image:: ../images/FAQBB.*
+``PropagatingModes`` contains the wavefunctions, velocities and momenta of the
+modes at the requested energy (2.5 in this example).  In order to understand
+the order in which these quantities are returned it is often useful to look at
+the a section of the band structure for the system in question:
 
-Here we take this Fermi energy so we can have 2 modes. Note that the velocity is negative for the incoming modes.
-We take the modes coming from the left lead (lead 0 in this case).
+.. image:: ../images/FAQPM.*
 
-.. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQAG
-    :end-before: #HIDDEN_END_FAQAG
+On the above band structure we have labelled the 4 modes in the order
+that they appear in the output of ``modes()`` at energy 2.5. Note that
+the modes are sorted in the following way:
 
-Then, we need to distinguish the different modes, it is indexed according to the decreasing k . For example here, we take the first mode and we want to know the wave function on the site at the position (6,2). To do that, we need to find the tag of this site using ``closest`` and then use ``id_by_site`` to find the index in the low-level system.
+    + First all the modes with negative velocity, then all the modes with
+      positive velocity
+    + Negative velocity modes are ordered by *increasing* momentum
+    + Positive velocity modes are ordered by *decreasing* momentum
 
-.. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQAH
-    :end-before: #HIDDEN_END_FAQAH
-
-
-How are degrees of freedom ordered ?
-======================================
+For more complicated systems and band structures this can lead to some
+possibly unintuitive orderings:
 
-We take the example from :doc:`here <../tutorial/superconductors>` which has 2 degrees of freedom. We will focus on the general ordering of the  degrees of freedom with the example of the ``wave_function()``
+.. image:: ../images/FAQPMC.*
 
-.. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQAO
-    :end-before: #HIDDEN_END_FAQAO
 
-.. image:: ../images/FAQTT.*
+How does Kwant order scattering states?
+=======================================
+Scattering states calculated using `~kwant.solvers.default.wave_function` are returned in the
+same order as the "incoming" modes of `~kwant.system.InfiniteSystem.modes`.
+Kwant considers that the translational symmetry of a lead points "towards
+infinity" (*not* towards the system) which means that the incoming modes are
+those that have *negative* velocities:
 
-.. image:: ../images/FAQSS.*
 
+How does Kwant order components of an individual wavefunction?
+==============================================================
+In `How can I get access to a system's sites?`_ we saw that the sites of a finalized system
+are available as a list through the ``sites`` attribute, and that one can look up the index
+of a site with the ``id_by_site`` attribute.
 
-With one degree of freedom, we had just to take the output of ``id_by_site`` to print the wave_function on a specific site.
+When all the site families present in a system have only 1 degree of freedom
+per site (i.e.  the all the onsites are scalars) then the index into a
+wavefunction defined over the system is exactly the site index:
 
 .. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQAH
-    :end-before: #HIDDEN_END_FAQAH
+    :start-after: #HIDDEN_BEGIN_ORD1
+    :end-before: #HIDDEN_END_ORD1
+.. literalinclude:: ../images/FAQORD1.txt
+
+We see that the wavefunction on a single site is a single complex number, as
+expected.
 
-In the case of 2 degrees of freedom, it changes the way we have access to a given site at a given orbital.
+If a site family have more than 1 degree of freedom per site (e.g. spin or
+particle-hole) then Kwant places degrees of freedom on the same site adjacent
+to one another.  In the case where all site families in the system have the
+*same* number of degrees of freedom, we can then simply *reshape* the
+wavefunction into a matrix, where the row number indexes the site, and the
+column number indexes the degree of freedom on that site:
 
 .. literalinclude:: FAQ.py
-    :start-after: #HIDDEN_BEGIN_FAQAP
-    :end-before: #HIDDEN_END_FAQAP
+    :start-after: #HIDDEN_BEGIN_ORD2
+    :end-before: #HIDDEN_END_ORD2
+.. literalinclude:: ../images/FAQORD2.txt
+
+We see that the wavefunction on a single site is a *vector* of 2 complex numbers,
+as we expect.
 
-The order of the orbitals and sites are completly modified. The N first orbitals of the first site are between the index 0 and N-1, for the second site, it is between N and 2N-1, etc...
+If there are different site families present in the system that have *different*
+numbers of orbitals per site, then the situation becomes much more involved,
+because we cannot simply "reshape" the wavefunction like we did in the
+preceding example.