diff --git a/doc/source/images/FAQ.py.diff b/doc/source/images/FAQ.py.diff
new file mode 100644
index 0000000000000000000000000000000000000000..b8649809d26ecbb1279ce6b781ce2bef4c13e45c
--- /dev/null
+++ b/doc/source/images/FAQ.py.diff
@@ -0,0 +1,223 @@
+--- original
++++ modified
+@@ -1,4 +1,4 @@
+-
++import _defs
+ from cmath import exp
+ import numpy as np
+ import kwant
+@@ -11,6 +11,13 @@
+
+ ######## 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)),
++ dpi=_defs.dpi, bbox_inches='tight')
++
+ a = 1
+ lat = kwant.lattice.square(a)
+ syst = kwant.Builder()
+@@ -19,6 +26,7 @@
+ syst[lat(1,1)] = 4
+
+ kwant.plot(syst)
++save_figure('FAQ122')
+
+
+
+@@ -40,7 +48,14 @@
+ 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')
+
+
+ ########## What is a hopping? #######################
+@@ -54,6 +69,7 @@
+ syst[(lat(1,0) , lat(1,1))] = -1
+
+ kwant.plot(syst)
++save_figure('FAQ124')
+
+
+
+@@ -73,7 +89,7 @@
+
+ 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')
+
+ ## Deleting 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')
+
+
+
+@@ -145,7 +162,7 @@
+
+ ## Plotting the system
+ plot_system(syst)
+-
++save_figure('FAQ6B')
+
+
+
+@@ -162,7 +179,7 @@
+
+ syst[kwant.builder.HoppingKind((1,0), lat)] = -1
+ kwant.plot(syst)
+-
++save_figure('FAQ3')
+ ## Polyatomic lattice
+
+ lat = kwant.lattice.kagome()
+@@ -193,12 +210,13 @@
+
+ syst[kwant.builder.HoppingKind((0,1), b, b)] = -1 # equivalent to syst[kwant.builder.HoppingKind((0,1), 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
+ plot_system(syst)
+-
++save_figure('FAQ11')
+
+
+ ########## How to create the hoppings between adjacent sites? ################
+@@ -213,11 +231,13 @@
+
+ syst[lat.neighbors()] = -1 ## equivalent to lat.neighbors(1)
+ kwant.plot(syst)
++save_figure('FAQ4')
+
+ del syst[lat.neighbors()] ## deletes every hoppings previously created to add new one
+ syst[lat.neighbors(2)] = -1
+
+ kwant.plot(syst)
++save_figure('FAQ5')
+
+ ## Polyatomic lattice
+
+@@ -236,14 +256,17 @@
+ ## .neighbors()
+ syst[lat.neighbors()] = -1
+ plot_system(syst)
++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
+
+
+ syst[lat.neighbors(2)] = -1
+ plot_system(syst)
++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
+ 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
+ plot_system(syst)
++save_figure('FAQAC')
+ ## Creation of the top lead
+
+ lat_lead = kwant.lattice.square(a)
+@@ -320,6 +345,8 @@
+
+ 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')
+
+
+
+@@ -356,6 +384,7 @@
+ cuboid.fill(model, cuboid_shape, (0, 0, 0));
+
+ kwant.plot(cuboid);
++save_figure('FAQaaa')
+
+ # Build electrode (black).
+ def electrode_shape(site):
+@@ -369,7 +398,7 @@
+
+ cuboid.attach_lead(electrode)
+ kwant.plot(cuboid);
+-
++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()
+
+@@ -401,6 +431,7 @@
+ plt.xlabel("momentum [(lattice constant)^-1]")
+ plt.ylabel("energy [t]")
+ plt.show()
++save_figure('FAQBB')
+ plt.clf()
+
+ wf = kwant.wave_function(fsyst, Ef)
+@@ -474,7 +505,7 @@
+
+ # 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')
+ plt.clf()
+
+ wf = kwant.wave_function(syst, Ef)
diff --git a/doc/source/tutorial/FAQ.py b/doc/source/tutorial/FAQ.py
new file mode 100644
index 0000000000000000000000000000000000000000..0df1b6a9c2e18b413daa69b639b0f312c1387d4d
--- /dev/null
+++ b/doc/source/tutorial/FAQ.py
@@ -0,0 +1,572 @@
+
+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 tinyarray
+
+
+######## What is a Site? ##############
+
+#HIDDEN_BEGIN_FAQ122
+a = 1
+lat = kwant.lattice.square(a)
+syst = kwant.Builder()
+
+syst[lat(1,0)] = 4
+syst[lat(1,1)] = 4
+
+kwant.plot(syst)
+#HIDDEN_END_FAQ122
+
+
+
+################# 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))
+
+## 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
+
+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
+#HIDDEN_BEGIN_FAQ124
+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:
+a = 2
+lat = kwant.lattice.cubic(a)
+syst = kwant.Builder()
+
+L = 10
+W = 10
+H = 2
+
+## Adding sites to the system:
+
+syst[ (lat(i,j,k) for i in range(L) for j in range(W) for k in range(H)) ] = 4
+kwant.plot(syst)
+
+
+## Deleting sites to create a hole
+
+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
+#HIDDEN_BEGIN_FAQ3
+## Before finalizing the system
+
+sites = []
+for site in builder.sites():
+ sites.append(site) ## here we choose to add the sites to a list
+
+#HIDDEN_END_FAQ3
+#HIDDEN_BEGIN_FAQ7
+## 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)
+#HIDDEN_END_FAQ7
+
+
+
+################ How to plot a polyatomic lattice with different colors? ##############"
+
+#HIDDEN_BEGIN_FAQ8
+lat = kwant.lattice.kagome()
+syst = kwant.Builder()
+
+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
+
+def plot_system(syst):
+
+ def family_color(site):
+ if site.family == a:
+ return 'red'
+ if site.family == b:
+ return 'green'
+ else:
+ return 'blue'
+
+ def hopping_lw(site1, site2):
+ return 0.1 if site1.family == site2.family else 0.1
+
+ 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[lat.neighbors()] = -1
+
+## Plotting the system
+plot_system(syst)
+
+#HIDDEN_END_FAQ9
+
+
+
+############### How to create every hoppings in a given direction using Hoppingkind? ################
+
+## Monatomic lattice
+
+#HIDDEN_BEGIN_FAQ4
+
+## Creation of 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[kwant.builder.HoppingKind((1,0), lat)] = -1
+kwant.plot(syst)
+#HIDDEN_END_FAQ4
+
+## Polyatomic lattice
+
+lat = kwant.lattice.kagome()
+syst = kwant.Builder()
+
+a , b, c = lat.sublattices ## The kagome lattice has 3 sublattices
+
+def plot_system(syst):
+
+ def family_color(site):
+ if site.family == a:
+ return 'blue'
+ if site.family == b:
+ return 'red'
+ else:
+ return 'green'
+
+
+ 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
+
+
+
+#HIDDEN_BEGIN_FAQ13
+syst[kwant.builder.HoppingKind((0,1), b, b)] = -1 # equivalent to syst[kwant.builder.HoppingKind((0,1), b)] = -1
+#HIDDEN_END_FAQ13
+plot_system(syst)
+del syst[kwant.builder.HoppingKind((0,1), b, b)] ## delete the hoppings previously created
+#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
+#HIDDEN_END_FAQ14
+plot_system(syst)
+
+
+
+########## How to create the hoppings between adjacent sites? ################
+
+## Monatomic lattice
+
+#HIDDEN_BEGIN_FAQ5
+
+## Creation of 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.neighbors()] = -1 ## equivalent to lat.neighbors(1)
+kwant.plot(syst)
+
+del syst[lat.neighbors()] ## deletes every hoppings previously created to add new one
+syst[lat.neighbors(2)] = -1
+
+kwant.plot(syst)
+#HIDDEN_END_FAQ5
+
+## Polyatomic lattice
+
+#HIDDEN_BEGIN_FAQ6
+
+## Hoppings using .neighbors()
+#HIDDEN_BEGIN_FAQ10
+## Creation of the system
+lat = kwant.lattice.kagome()
+syst = kwant.Builder()
+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
+
+## .neighbors()
+syst[lat.neighbors()] = -1
+#HIDDEN_END_FAQ10
+plot_system(syst)
+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
+
+
+#HIDDEN_BEGIN_FAQ12A
+syst[lat.neighbors(2)] = -1
+#HIDDEN_END_FAQ12A
+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
+
+def plot_system(syst):
+
+ def family_color(site):
+ if site.family == subA:
+ return 'blue'
+ if site.family == subB:
+ return 'yellow'
+ 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()
+
+L = 5
+W = 5
+
+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[lat.neighbors()] = -1
+plot_system(syst)
+#HIDDEN_END_FAQAB
+#HIDDEN_BEGIN_FAQAC
+## We manually add sites of the same lead lattice
+
+
+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)
+
+
+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
+
+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.
+cubic = kwant.lattice.cubic()
+sym_3d = kwant.TranslationalSymmetry([1, 0, 0], [0, 1, 0], [0, 0, 1])
+model = kwant.Builder(sym_3d)
+model[cubic(0, 0, 0)] = 4
+model[cubic.neighbors()] = -1
+#HIDDEN_END_FAQccc
+
+#HIDDEN_BEGIN_FAQddd
+# Build scattering region (white).
+def cuboid_shape(site):
+ x, y, z = abs(site.pos)
+ return x < 4 and y < 10 and z < 3
+
+cuboid = kwant.Builder()
+cuboid.fill(model, cuboid_shape, (0, 0, 0));
+
+kwant.plot(cuboid);
+#HIDDEN_END_FAQddd
+
+#HIDDEN_BEGIN_FAQeee
+# Build electrode (black).
+def electrode_shape(site):
+ x, y, z = site.pos - (0, 5, 2)
+ 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
+
+cuboid.fill(electrode, lambda s: abs(s.pos[0]) < 7, (0, 5, 4)); ## scattering region
+
+cuboid.attach_lead(electrode)
+kwant.plot(cuboid);
+#HIDDEN_END_FAQeee
+
+
+###### How to extract the wavefunction informations on a specific site? ###############
+
+#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)
+
+sym_lead = kwant.TranslationalSymmetry((-a,0))
+lead = kwant.Builder(sym_lead)
+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
+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
+
+
+#HIDDEN_BEGIN_FAQAH
+
+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 = 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
+
+
+#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)
+
+
+ ### 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.
+ syst = kwant.Builder()
+
+ #### 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()
+
+wf = kwant.wave_function(syst, Ef)
+wf_left_lead = wf(0) # Wave function for the first lead (left)
+
+
+
+#HIDDEN_END_FAQAO
+
+
+#HIDDEN_BEGIN_FAQAP
+nb_degrees = 2 ## 2 degrees of freedom
+
+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 = []
+
+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
+
+
+#HIDDEN_END_FAQAP
diff --git a/doc/source/tutorial/FAQ.rst b/doc/source/tutorial/FAQ.rst
new file mode 100644
index 0000000000000000000000000000000000000000..48ec9abadb011ca18846b6c0e0332a178e3f9614
--- /dev/null
+++ b/doc/source/tutorial/FAQ.rst
@@ -0,0 +1,335 @@
+---------------------------------------------------------------
+Frequently Asked Questions :
+---------------------------------------------------------------
+
+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 site?
+=================
+
+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 :
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ122
+ :end-before: #HIDDEN_END_FAQ122
+
+.. 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``
+
+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
+
+
+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.
+
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ123
+ :end-before: #HIDDEN_END_FAQ123
+
+.. 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.
+
+What is a family ? a tag?
+============================
+
+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.
+
+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 :
+
+.. 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.
+
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ3
+ :end-before: #HIDDEN_END_FAQ3
+
+Before finalizing, we get access to the sites with : ``syst.sites()`` .
+
+.. 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.
+
+How to plot a polyatomic lattice with different colors?
+===============================================================================
+We take a kagome lattice as example here, it has 3 sublattices.
+
+.. 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.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ9
+ :end-before: #HIDDEN_END_FAQ9
+
+.. image:: ../images/FAQ6B.*
+
+How to create every hoppings in a given direction using ``Hoppingkind``?
+==========================================================================
+
+
+If the lattice is monatomic:
+
+``HoppingKind`` can be used to easily create a hopping in a specific direction.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ4
+ :end-before: #HIDDEN_END_FAQ4
+
+.. 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.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ13
+ :end-before: #HIDDEN_END_FAQ13
+
+Here, we want the hoppings between the sites from sublattice b with a direction of (0,1) in the lattice coordinates.
+
+.. image:: ../images/FAQ10.*
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ14
+ :end-before: #HIDDEN_END_FAQ14
+
+Here, we create hoppings between the sites of the same lattice coordinates but from different families.
+
+.. image:: ../images/FAQ11.*
+
+
+
+How to 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.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ5
+ :end-before: #HIDDEN_END_FAQ5
+
+.. image:: ../images/FAQ4.*
+.. image:: ../images/FAQ5.*
+
+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.
+
+.. 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.*
+
+.. 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.*
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQ12A
+ :end-before: #HIDDEN_END_FAQ12A
+
+It creates hoppings between the second nearest neighbors of every sites of the system.
+
+.. image:: ../images/FAQ9.*
+
+
+
+How to create a lead with a lattice different from the scattering region?
+===========================================================================
+
+First of all, we need to plot the sites in different colors depending on the sites families.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAA
+ :end-before: #HIDDEN_END_FAQAA
+
+
+Then, we create the scattering region, here we take the example of the honeycomb lattice.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAB
+ :end-before: #HIDDEN_END_FAQAB
+
+.. 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.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAC
+ :end-before: #HIDDEN_END_FAQAC
+
+.. 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.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAD
+ :end-before: #HIDDEN_END_FAQAD
+
+.. 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.
+
+
+We first create the system template with its symmetries.
+
+.. 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
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQddd
+ :end-before: #HIDDEN_END_FAQddd
+
+.. image:: ../images/FAQaaa.*
+
+Finally, we create the template that will used to create the second part of the scattering region and the lead.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQeee
+ :end-before: #HIDDEN_END_FAQeee
+
+.. 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.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAF
+ :end-before: #HIDDEN_END_FAQAF
+
+.. image:: ../images/FAQBA.*
+.. image:: ../images/FAQBB.*
+
+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).
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAG
+ :end-before: #HIDDEN_END_FAQAG
+
+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.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAH
+ :end-before: #HIDDEN_END_FAQAH
+
+
+How are degrees of freedom ordered ?
+======================================
+
+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()``
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAO
+ :end-before: #HIDDEN_END_FAQAO
+
+.. image:: ../images/FAQTT.*
+
+.. image:: ../images/FAQSS.*
+
+
+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.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAH
+ :end-before: #HIDDEN_END_FAQAH
+
+In the case of 2 degrees of freedom, it changes the way we have access to a given site at a given orbital.
+
+.. literalinclude:: FAQ.py
+ :start-after: #HIDDEN_BEGIN_FAQAP
+ :end-before: #HIDDEN_END_FAQAP
+
+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...
diff --git a/doc/source/tutorial/index.rst b/doc/source/tutorial/index.rst
index 2cade40aa1eec20072abad7ae1360141af1e8569..bee5d10addfc6969f7164b0b2b3bfc56e5c77296 100644
--- a/doc/source/tutorial/index.rst
+++ b/doc/source/tutorial/index.rst
@@ -12,3 +12,4 @@ Tutorial: learning Kwant through examples
plotting
kpm
discretize
+ FAQ