diff --git a/doc/source/tutorial/faq.rst b/doc/source/tutorial/faq.rst
index 3833f4743bdd9450f051a426c3ad745ee9005fa1..1022bf03982fa419a95faba812013ee02f2771d4 100644
--- a/doc/source/tutorial/faq.rst
+++ b/doc/source/tutorial/faq.rst
@@ -434,7 +434,12 @@ Scattering states calculated using `~kwant.solvers.default.wave_function` are re
 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:
+those that have *negative* velocities.
+
+This means that for a lead attached on the left of a scattering region (with
+symmetry vector $(-1, 0)$, for example), the
+positive $k$ direction (when inspecting the lead's band structure) actually
+corresponds to the *negative* $x$ direction.
 
 
 How does Kwant order components of an individual wavefunction?
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index 61883c43efafc8dbca663f415bdc9c96d99d1195..222ad0ebb84400ceaa7dc7cfa2ec705f5ef3e955 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -135,6 +135,9 @@ class PropagatingModes:
     velocity, `k` is momentum and `conserved_quantity` is the conservation
     law eigenvalue.
 
+    In the above, the positive velocity and momentum directions are defined
+    with respect to the translational symmetry direction of the system.
+
     The first dimension of `wave_functions` corresponds to the orbitals of all
     the sites in a unit cell, the second one to the number of the mode.  Each
     mode is normalized to carry unit current. If several modes have the same
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index f5f6197e9b791b2fb02c236ed0e26153fdc93f97..fbe88b69219d71e4a68f828bf306b5ef2381864a 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -570,13 +570,23 @@ class SparseSolver(metaclass=abc.ABCMeta):
 
         Notes
         -----
-
         The returned object can be itself called like a function.  Given a lead
         number, it returns a 2d NumPy array that contains the wave function
         within the scattering region due to each incoming mode of the given
-        lead.  Index 0 is the mode number, index 1 is the orbital number.  The
-        modes appear in the same order as incoming modes in
-        `kwant.physics.modes`.
+        lead.  Index 0 is the mode number, index 1 is the orbital number.
+
+        The modes appear in the same order as the negative velocity modes in
+        `kwant.physics.PropagatingModes`. In Kwant's convention leads are attached
+        so that their translational symmetry points *away* from the scattering
+        region::
+
+             left lead    SR   right lead
+             /---------\ /---\ /---------\
+             ...-3-2-1-0-X-X-X-0-1-2-3-...
+
+        This means that incoming modes (coming from infinity towards the
+        scattering region) have *negative* velocity with respect to the
+        lead's symmetry direction.
 
         Examples
         --------
diff --git a/kwant/system.py b/kwant/system.py
index 62158d41819ef8f2b62a3e2e8c69f1ceed10f0f6..dfa28558d0efdb08a7ef8af788299e474e6b2759 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -224,6 +224,12 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
 
         See documentation of `~kwant.physics.PropagatingModes` and
         `~kwant.physics.StabilizedModes` for the return format details.
+
+        The wave functions of the returned modes are defined over the
+        *unit cell* of the system, which corresponds to the degrees of
+        freedom on the first ``cell_sites`` sites of the system
+        (recall that infinite systems store first the sites in the unit
+        cell, then connected sites in the neighboring unit cell).
         """
         from . import physics   # Putting this here avoids a circular import.
         ham = self.cell_hamiltonian(args, params=params)