diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index f52d5c630b84008acf856f41f82cdf99b4218028..ee2257104f50e9f0ecad80210b108039c31bb4e2 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -57,18 +57,27 @@ check for dependencies installed:
 build documentation:
   stage: test
   script:
-    - pip3 install sympy
     - make -C doc realclean; make -C doc html SPHINXOPTS='-A website_deploy=True -n -W' REFNAME="${CI_COMMIT_TAG:-$CI_COMMIT_SHA}" SOURCE_URL="$CI_PROJECT_URL"/blob
   artifacts:
     paths:
       - doc/build/html/
     expire_in: 1 month
 
+build PDF documentation:
+  stage: test
+  script:
+    - pip3 install sympy
+    - make -C doc latex SPHINXOPTS='-n -W'
+    - cd doc/build/latex
+    - make all-pdf
+  artifacts:
+    paths:
+      - doc/build/latex/kwant.pdf
+    expire_in: 1 month
 
 run tests:
   stage: test
   script:
-    - pip3 install sympy
     - py.test -r w --cov=kwant --cov-report term --cov-report html --flakes kwant
   artifacts:
     paths:
@@ -78,7 +87,6 @@ run tests:
 check for broken links in doc:
   stage: test
   script:
-    - pip3 install sympy
     - make -C doc linkcheck
   allow_failure: true
 
diff --git a/doc/Makefile b/doc/Makefile
index b9fd88395b2ca3a84cd4de24648bbd49ff8914c7..1da6042002fa775ca6dbd688a65a8d6bfdf95b4f 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -127,7 +127,7 @@ doctest:
 	      "results in $(BUILDDIR)/doctest/output.txt."
 
 %.pdf: %.svg
-	rsvg-convert -f pdf -o $@ $<
+	inkscape --export-pdf=$@ $<
 
 # Make tutorial scripts by extracting the (complete!) context of the "patches".
 # We make sure not to use 'wiggle' here.
diff --git a/doc/source/code/figure/magnetic_texture.py.diff b/doc/source/code/figure/magnetic_texture.py.diff
index 0b1f90cb6bd17d88a97b7e1924641fb640b15a63..a29a458e43ea5314b09e93426bf907872228f8f9 100644
--- a/doc/source/code/figure/magnetic_texture.py.diff
+++ b/doc/source/code/figure/magnetic_texture.py.diff
@@ -40,7 +40,7 @@
      theta = (tanh(r_tilde) - 1) * (pi / 2)
  
      if r == 0:
-         m_i = [0, 0, 1]
+         m_i = [0, 0, -1]
      else:
          m_i = [
              (x / r) * sin(theta),
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 9b5c14cd31efe65685d119bee5b305def1a160f4..f0272556a1cee688b73c60fdea6574d8792ba5ee 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -199,6 +199,8 @@ r"""\makeatletter
 \newcommand{\braket}[2]{\left\langle#1|#2\right\rangle}
 \newcommand{\ri}{\text{i}}
 \newcommand{\rd}{\text{d}}
+
+\usepackage{unicode-math}
 """}
 
 # Grouping the document tree into LaTeX files. List of tuples
@@ -210,6 +212,8 @@ latex_documents = [
    'manual'),
 ]
 
+latex_engine = 'xelatex'
+
 # The name of an image file (relative to this directory) to place at the top of
 # the title page.
 #latex_logo = None
diff --git a/doc/source/pre/whatsnew/1.3.rst b/doc/source/pre/whatsnew/1.3.rst
index 0f4e250fdcf608271da94544e6396cddd4638ea8..707608cb65e2d54c7854f01b0ded189c2ce2d815 100644
--- a/doc/source/pre/whatsnew/1.3.rst
+++ b/doc/source/pre/whatsnew/1.3.rst
@@ -52,8 +52,8 @@ the following::
     psi = kwant.wave_function(syst)(0)[0]
 
     # create the operators
-    Q = kwant.physics.LocalOperator(syst)
-    J = kwant.physics.Current(syst)
+    Q = kwant.operator.Density(syst)
+    J = kwant.operator.Current(syst)
 
     # evaluate the expectation value with the wavefunction
     q = Q(psi)
diff --git a/doc/source/tutorial/discretize.rst b/doc/source/tutorial/discretize.rst
index c4a5b6a864db1843470365510ee8c73d44fe5d78..dd8f65386b2477a837dd7663a6b9c23bd31160d1 100644
--- a/doc/source/tutorial/discretize.rst
+++ b/doc/source/tutorial/discretize.rst
@@ -219,7 +219,7 @@ and its discretized approximation
 
 
 where :math:`t=\frac{\hbar^2}{2ma^2}`, are only valid in the limit
-:math:`E \lt t`. The grid spacing :math:`a` must be chosen according
+:math:`E < t`. The grid spacing :math:`a` must be chosen according
 to how high in energy you need your tight-binding model to be valid.
 
 It is possible to set :math:`a` through the ``grid_spacing`` parameter
diff --git a/doc/source/tutorial/operators.rst b/doc/source/tutorial/operators.rst
index c7e66feaf40efcd00339571bef2efa462d217838..99fef17fe0deeda4d9fe7a5b71e17c4885e0dcc8 100644
--- a/doc/source/tutorial/operators.rst
+++ b/doc/source/tutorial/operators.rst
@@ -34,12 +34,12 @@ that, on site :math:`i`, points in the direction of the unit vector
 We shall take the following form for :math:`\mathbf{m}_i`:
 
 .. math::
-    \mathbf{m}_i &=&\ \left(
+    \mathbf{m}_i &=\ \left(
         \frac{x_i}{x_i^2 + y_i^2} \sin θ_i,\
         \frac{y_i}{x_i^2 + y_i^2} \sin θ_i,\
         \cos θ_i \right)^T,
     \\
-    θ_i &=&\ \frac{π}{2} \tanh \frac{r_i - r_0}{δ},
+    θ_i &=\ \frac{π}{2} (\tanh \frac{r_i - r_0}{δ} - 1),
 
 where :math:`x_i` and :math:`y_i` are the :math:`x` and :math:`y` coordinates
 of site :math:`i`, and :math:`r_i = \sqrt{x_i^2 + y_i^2}`.
diff --git a/kwant/builder.py b/kwant/builder.py
index 997a7e9cbaa2810c574b1415f90100950324e139..5a6737c2b73a1874fc9ce39373df3dd103f9fb72 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -1426,7 +1426,10 @@ class Builder:
     def attach_lead(self, lead_builder, origin=None, add_cells=0):
         """Attach a lead to the builder, possibly adding missing sites.
 
-        Returns the lead number (integer) of the attached lead.
+        This method first adds sites from 'lead_builder' until the interface
+        where the lead will attach is "smooth".  Then it appends the
+        'lead_builder' and the interface sites as a
+        `~kwant.builder.BuilderLead` to the 'leads' of this builder.
 
         Parameters
         ----------
@@ -1453,11 +1456,41 @@ class Builder:
 
         Notes
         -----
-        This method is not fool-proof, i.e. if it returns an error, there is
+        This method is not fool-proof, i.e. if it raises an error, there is
         no guarantee that the system stayed unaltered.
 
-        The lead numbering starts from zero and increments from there, i.e.
-        the leads are numbered in the order in which they are attached.
+        The system must "interrupt" the lead that is being attached. This means
+        that for each site in the lead that has a hopping to a neighbnoring
+        unit cell there must be at least one site in the system that is an
+        image of the lead site under the action of the lead's translational
+        symmetry.  In order to interrupt the lead, the system must contain
+        sites from the same site family as the sites in the lead.  Below are
+        three examples of leads being attached to systems::
+
+            Successful           Successful          Unsuccessful
+            ----------           ----------          ------------
+            Lead   System        Lead   System        Lead  System
+            x-  …      o         x   …                x-  …
+            |          |         |                    |
+            x-  …    o-o         x-  …    o-o         x-  …  o-o
+            |        | |         |        | |         |      | |
+            x-  …  o-o-o         x-  …  o-o-o         x-  …  o-o
+
+        The second case succeeds, as even though the top site has no image in
+        the system, because the top site has no hoppings to sites in other unit
+        cells.
+
+        Sites may be added to the system when the lead is attached, so that the
+        interface to the lead is "smooth". Below we show the system after
+        having attached a lead. The 'x' symbols in the system indicate the
+        added sites::
+
+            Lead   System        Lead   System
+            x-  …  x-x-o         x   …  x
+            |      | | |         |      |
+            x-  …  x-o-o         x-  …  x-o-o
+            |      | | |         |      | | |
+            x-  …  o-o-o         x-  …  o-o-o
         """
         if self.symmetry.num_directions:
             raise ValueError("Leads can only be attached to finite systems.")
diff --git a/kwant/graph/core.pyx b/kwant/graph/core.pyx
index 1f9c15f2249d7a186eb2134e3bed4b0d9d930dcc..a5175bee100b26328fc4fa8fd516799c5f0fa14b 100644
--- a/kwant/graph/core.pyx
+++ b/kwant/graph/core.pyx
@@ -719,15 +719,22 @@ cdef class CGraph_malloc(CGraph):
         init_args = (twoway, edge_nr_translation, num_nodes,
                 num_pp_edges, num_pn_edges, num_np_edges)
 
-        return (init_args, self._heads_idxs, self._heads, self._tails_idxs,
-                self._tails, self._edge_ids, self._edge_ids_by_edge_nr)
+        return init_args, (self._heads_idxs, self._heads, self._tails_idxs,
+                           self._tails, self._edge_ids,
+                           self._edge_ids_by_edge_nr)
 
     def __setstate__(self, state):
-        self.__init__(*state[0])
+        init_args, arrays = state
+        self.__init__(*init_args)
         array_attributes = (self._heads_idxs, self._heads, self._tails_idxs,
                             self._tails, self._edge_ids,
                             self._edge_ids_by_edge_nr)
-        for attribute, value in zip(array_attributes, state[1:]):
+        for attribute, value in zip(array_attributes, arrays):
             if attribute is None:
                 continue
             attribute[:] = value
+
+    # We are required to implement this as of Cython 0.26
+    def __reduce__(self):
+        state = init_args, _ = self.__getstate__()
+        return (CGraph_malloc, init_args, state, None, None)
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index f2ac1a8ebc3f20a96e7aa4b043a24a3067c8a222..4fcdd73ba925a8ff27cf0071cd4fc48ee4d4728c 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -116,6 +116,10 @@ class PropagatingModes:
     eigenvalue the modes with negative velocity are ordered by increasing
     momentum, and the modes with positive velocity are ordered by decreasing
     momentum. Finally, modes are ordered by the magnitude of their velocity.
+    To summarize, the modes are ordered according to the key
+    `(sign(v), conserved_quantity, sign(v) * k , abs(v))` where `v` is
+    velocity, `k` is momentum and `conserved_quantity` is the conservation
+    law eigenvalue.
 
     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
@@ -1013,12 +1017,11 @@ def modes(h_cell, h_hop, tol=1e6, stabilization=None, *,
 
     Notes
     -----
-    The propagating modes are sorted according to the longitudinal component of
-    their k-vector, with incoming modes having k sorted in descending order,
-    and outgoing modes having k sorted in ascending order.  In simple cases
-    where bands do not cross, this ordering corresponds to "lowest modes
-    first". In general, however, it is necessary to examine the band structure
-    -- something this function is not doing by design.
+    The sorting of the propagating modes is fully described in the
+    documentation for `~kwant.physics.PropagatingModes`.  In simple cases where
+    bands do not cross, this ordering corresponds to "lowest modes first". In
+    general, however, it is necessary to examine the band structure --
+    something this function is not doing by design.
 
     Propagating modes with the same momentum are orthogonalized. All the
     propagating modes are normalized by current.