diff --git a/doc/source/images/magnetic_texture.py.diff b/doc/source/images/magnetic_texture.py.diff
new file mode 100644
index 0000000000000000000000000000000000000000..bc7acb4dcf1b7bda417c3dc3458eb61655435e05
--- /dev/null
+++ b/doc/source/images/magnetic_texture.py.diff
@@ -0,0 +1,143 @@
+--- original
++++ modified
+@@ -11,6 +11,9 @@
+ # - operators
+ # - plotting vector fields
+
++import _defs
++from contextlib import redirect_stdout
++
+ from math import sin, cos, tanh, pi
+ import itertools
+ import numpy as np
+@@ -79,7 +82,13 @@
+ return syst
+
+
+-def plot_vector_field(syst, params):
++def save_plot(fname):
++ for extension in ('.pdf', '.png'):
++ plt.savefig(fname + extension,
++ dpi=_defs.dpi, bbox_inches='tight')
++
++
++def plot_vector_field(syst, params, fname=None):
+ xmin, ymin = min(s.tag for s in syst.sites)
+ xmax, ymax = max(s.tag for s in syst.sites)
+ x, y = np.meshgrid(np.arange(xmin, xmax+1), np.arange(ymin, ymax+1))
+@@ -88,28 +97,38 @@
+ m_i = np.reshape(m_i, x.shape + (3,))
+ m_i = np.rollaxis(m_i, 2, 0)
+
+- fig, ax = plt.subplots(1, 1)
++ fig, ax = plt.subplots(1, 1, figsize=(9, 7))
+ im = ax.quiver(x, y, *m_i, pivot='mid', scale=75)
+ fig.colorbar(im)
+- plt.show()
++ if fname:
++ save_plot(fname)
++ else:
++ plt.show()
+
+
+-def plot_densities(syst, densities):
+- fig, axes = plt.subplots(1, len(densities))
++def plot_densities(syst, densities, fname=None):
++ fig, axes = plt.subplots(1, len(densities), figsize=(7*len(densities), 7))
+ for ax, (title, rho) in zip(axes, densities):
+ kwant.plotter.map(syst, rho, ax=ax, a=4)
+ ax.set_title(title)
+- plt.show()
++ if fname:
++ save_plot(fname)
++ else:
++ plt.show()
++
+
+
+-def plot_currents(syst, currents):
+- fig, axes = plt.subplots(1, len(currents))
++def plot_currents(syst, currents, fname=None):
++ fig, axes = plt.subplots(1, len(currents), figsize=(7*len(currents), 7))
+ if not hasattr(axes, '__len__'):
+ axes = (axes,)
+ for ax, (title, current) in zip(axes, currents):
+ kwant.plotter.current(syst, current, ax=ax, colorbar=False)
+ ax.set_title(title)
+- plt.show()
++ if fname:
++ save_plot(fname)
++ else:
++ plt.show()
+
+
+ def main():
+@@ -119,7 +138,7 @@
+ wf = kwant.wave_function(syst, energy=-1, params=params)
+ psi = wf(0)[0]
+
+- plot_vector_field(syst, dict(r0=20, delta=10))
++ plot_vector_field(syst, dict(r0=20, delta=10), fname='mag_field_direction')
+
+ # even (odd) indices correspond to spin up (down)
+ up, down = psi[::2], psi[1::2]
+@@ -144,7 +163,7 @@
+ ('$σ_0$', density),
+ ('$σ_z$', spin_z),
+ ('$σ_y$', spin_y),
+- ])
++ ], fname='spin_densities')
+
+ J_0 = kwant.operator.Current(syst)
+ J_z = kwant.operator.Current(syst, sigma_z)
+@@ -159,7 +178,7 @@
+ ('$J_{σ_0}$', current),
+ ('$J_{σ_z}$', spin_z_current),
+ ('$J_{σ_y}$', spin_y_current),
+- ])
++ ], fname='spin_currents')
+
+ def following_m_i(site, r0, delta):
+ m_i = field_direction(site.pos, r0, delta)
+@@ -173,7 +192,7 @@
+ plot_currents(syst, [
+ (r'$J_{\mathbf{m}_i}$', m_current),
+ ('$J_{σ_z}$', spin_z_current),
+- ])
++ ], fname='spin_current_comparison')
+
+
+ def circle(site):
+@@ -183,7 +202,9 @@
+
+ all_states = np.vstack((wf(0), wf(1)))
+ dos_in_circle = sum(rho_circle(p) for p in all_states) / (2 * pi)
+- print('density of states in circle:', dos_in_circle)
++ with open('circle_dos.txt', 'w') as f:
++ with redirect_stdout(f):
++ print('density of states in circle:', dos_in_circle)
+
+ def left_cut(site_to, site_from):
+ return site_from.pos[0] <= -39 and site_to.pos[0] > -39
+@@ -197,8 +218,10 @@
+ Jz_left = kwant.operator.Current(syst, sigma_z, where=left_cut, sum=True)
+ Jz_right = kwant.operator.Current(syst, sigma_z, where=right_cut, sum=True)
+
+- print('J_left:', J_left(psi), ' J_right:', J_right(psi))
+- print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
++ with open('current_cut.txt', 'w') as f:
++ with redirect_stdout(f):
++ print('J_left:', J_left(psi), ' J_right:', J_right(psi))
++ print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
+
+ J_m = kwant.operator.Current(syst, following_m_i)
+ J_z = kwant.operator.Current(syst, sigma_z)
+@@ -214,7 +237,7 @@
+ plot_currents(syst, [
+ (r'$J_{\mathbf{m}_i}$ (from left)', J_m_left),
+ (r'$J_{σ_z}$ (from left)', J_z_left),
+- ])
++ ], fname='bound_current')
+
+
+ if __name__ == '__main__':
diff --git a/doc/source/tutorial/index.rst b/doc/source/tutorial/index.rst
index 1729a23126f5c0a7cfbd46aaa94e44a64881adf6..275e660e1c10c58d806422e917b77b17896dad5a 100644
--- a/doc/source/tutorial/index.rst
+++ b/doc/source/tutorial/index.rst
@@ -11,3 +11,4 @@ Tutorial: learning Kwant through examples
tutorial6
tutorial7
tutorial8
+ tutorial9
diff --git a/doc/source/tutorial/magnetic_texture.py b/doc/source/tutorial/magnetic_texture.py
new file mode 100644
index 0000000000000000000000000000000000000000..89c7acfb793aecce1126a50fa45f31b360c2e053
--- /dev/null
+++ b/doc/source/tutorial/magnetic_texture.py
@@ -0,0 +1,245 @@
+# Tutorial 2.7. Spin textures
+# ===========================
+#
+# Physics background
+# ------------------
+# - Spin textures
+# - Skyrmions
+#
+# Kwant features highlighted
+# --------------------------
+# - operators
+# - plotting vector fields
+
+from math import sin, cos, tanh, pi
+import itertools
+import numpy as np
+import tinyarray as ta
+import matplotlib.pyplot as plt
+
+import kwant
+
+sigma_0 = ta.array([[1, 0], [0, 1]])
+sigma_x = ta.array([[0, 1], [1, 0]])
+sigma_y = ta.array([[0, -1j], [1j, 0]])
+sigma_z = ta.array([[1, 0], [0, -1]])
+
+# vector of Pauli matrices σ_αiβ where greek
+# letters denote spinor indices
+sigma = np.rollaxis(np.array([sigma_x, sigma_y, sigma_z]), 1)
+
+#HIDDEN_BEGIN_model
+def field_direction(pos, r0, delta):
+ x, y = pos
+ r = np.linalg.norm(pos)
+ r_tilde = (r - r0) / delta
+ theta = (tanh(r_tilde) - 1) * (pi / 2)
+
+ if r == 0:
+ m_i = [0, 0, 1]
+ else:
+ m_i = [
+ (x / r) * sin(theta),
+ (y / r) * sin(theta),
+ cos(theta),
+ ]
+
+ return np.array(m_i)
+
+
+def scattering_onsite(site, r0, delta, J):
+ m_i = field_direction(site.pos, r0, delta)
+ return J * np.dot(m_i, sigma)
+
+
+def lead_onsite(site, J):
+ return J * sigma_z
+#HIDDEN_END_model
+
+
+#HIDDEN_BEGIN_syst
+lat = kwant.lattice.square(norbs=2)
+
+def make_system(L=80):
+
+ syst = kwant.Builder()
+
+ def square(pos):
+ return all(-L/2 < p < L/2 for p in pos)
+
+ syst[lat.shape(square, (0, 0))] = scattering_onsite
+ syst[lat.neighbors()] = -sigma_0
+
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)),
+ conservation_law=-sigma_z)
+
+ lead[lat.shape(square, (0, 0))] = lead_onsite
+ lead[lat.neighbors()] = -sigma_0
+
+ syst.attach_lead(lead)
+ syst.attach_lead(lead.reversed())
+
+ return syst
+#HIDDEN_END_syst
+
+
+def plot_vector_field(syst, params):
+ xmin, ymin = min(s.tag for s in syst.sites)
+ xmax, ymax = max(s.tag for s in syst.sites)
+ x, y = np.meshgrid(np.arange(xmin, xmax+1), np.arange(ymin, ymax+1))
+
+ m_i = [field_direction(p, **params) for p in zip(x.flat, y.flat)]
+ m_i = np.reshape(m_i, x.shape + (3,))
+ m_i = np.rollaxis(m_i, 2, 0)
+
+ fig, ax = plt.subplots(1, 1)
+ im = ax.quiver(x, y, *m_i, pivot='mid', scale=75)
+ fig.colorbar(im)
+ plt.show()
+
+
+def plot_densities(syst, densities):
+ fig, axes = plt.subplots(1, len(densities))
+ for ax, (title, rho) in zip(axes, densities):
+ kwant.plotter.map(syst, rho, ax=ax, a=4)
+ ax.set_title(title)
+ plt.show()
+
+
+def plot_currents(syst, currents):
+ fig, axes = plt.subplots(1, len(currents))
+ if not hasattr(axes, '__len__'):
+ axes = (axes,)
+ for ax, (title, current) in zip(axes, currents):
+ kwant.plotter.current(syst, current, ax=ax, colorbar=False)
+ ax.set_title(title)
+ plt.show()
+
+
+def main():
+ syst = make_system().finalized()
+
+#HIDDEN_BEGIN_wavefunction
+ params = dict(r0=20, delta=10, J=1)
+ wf = kwant.wave_function(syst, energy=-1, params=params)
+ psi = wf(0)[0]
+#HIDDEN_END_wavefunction
+
+ plot_vector_field(syst, dict(r0=20, delta=10))
+
+#HIDDEN_BEGIN_ldos
+ # even (odd) indices correspond to spin up (down)
+ up, down = psi[::2], psi[1::2]
+ density = np.abs(up)**2 + np.abs(down)**2
+#HIDDEN_END_ldos
+
+#HIDDEN_BEGIN_lsdz
+ # spin down components have a minus sign
+ spin_z = np.abs(up)**2 - np.abs(down)**2
+#HIDDEN_END_lsdz
+
+#HIDDEN_BEGIN_lsdy
+ # spin down components have a minus sign
+ spin_y = 1j * (down.conjugate() * up - up.conjugate() * down)
+#HIDDEN_END_lsdy
+
+#HIDDEN_BEGIN_lden
+ rho = kwant.operator.Density(syst)
+ rho_sz = kwant.operator.Density(syst, sigma_z)
+ rho_sy = kwant.operator.Density(syst, sigma_y)
+
+ # calculate the expectation values of the operators with 'psi'
+ density = rho(psi)
+ spin_z = rho_sz(psi)
+ spin_y = rho_sy(psi)
+#HIDDEN_END_lden
+
+ plot_densities(syst, [
+ ('$σ_0$', density),
+ ('$σ_z$', spin_z),
+ ('$σ_y$', spin_y),
+ ])
+
+#HIDDEN_BEGIN_current
+ J_0 = kwant.operator.Current(syst)
+ J_z = kwant.operator.Current(syst, sigma_z)
+ J_y = kwant.operator.Current(syst, sigma_y)
+
+ # calculate the expectation values of the operators with 'psi'
+ current = J_0(psi)
+ spin_z_current = J_z(psi)
+ spin_y_current = J_y(psi)
+#HIDDEN_END_current
+
+ plot_currents(syst, [
+ ('$J_{σ_0}$', current),
+ ('$J_{σ_z}$', spin_z_current),
+ ('$J_{σ_y}$', spin_y_current),
+ ])
+
+#HIDDEN_BEGIN_following
+ def following_m_i(site, r0, delta):
+ m_i = field_direction(site.pos, r0, delta)
+ return np.dot(m_i, sigma)
+
+ J_m = kwant.operator.Current(syst, following_m_i)
+
+ # evaluate the operator
+ m_current = J_m(psi, params=dict(r0=25, delta=10))
+#HIDDEN_END_following
+
+ plot_currents(syst, [
+ (r'$J_{\mathbf{m}_i}$', m_current),
+ ('$J_{σ_z}$', spin_z_current),
+ ])
+
+
+#HIDDEN_BEGIN_density_cut
+ def circle(site):
+ return np.linalg.norm(site.pos) < 20
+
+ rho_circle = kwant.operator.Density(syst, where=circle, sum=True)
+
+ all_states = np.vstack((wf(0), wf(1)))
+ dos_in_circle = sum(rho_circle(p) for p in all_states) / (2 * pi)
+ print('density of states in circle:', dos_in_circle)
+#HIDDEN_END_density_cut
+
+#HIDDEN_BEGIN_current_cut
+ def left_cut(site_to, site_from):
+ return site_from.pos[0] <= -39 and site_to.pos[0] > -39
+
+ def right_cut(site_to, site_from):
+ return site_from.pos[0] < 39 and site_to.pos[0] >= 39
+
+ J_left = kwant.operator.Current(syst, where=left_cut, sum=True)
+ J_right = kwant.operator.Current(syst, where=right_cut, sum=True)
+
+ Jz_left = kwant.operator.Current(syst, sigma_z, where=left_cut, sum=True)
+ Jz_right = kwant.operator.Current(syst, sigma_z, where=right_cut, sum=True)
+
+ print('J_left:', J_left(psi), ' J_right:', J_right(psi))
+ print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
+#HIDDEN_END_current_cut
+
+#HIDDEN_BEGIN_bind
+ J_m = kwant.operator.Current(syst, following_m_i)
+ J_z = kwant.operator.Current(syst, sigma_z)
+
+ J_m_bound = J_m.bind(params=dict(r0=25, delta=10, J=1))
+ J_z_bound = J_z.bind(params=dict(r0=25, delta=10, J=1))
+
+ # Sum current local from all scattering states on the left at energy=-1
+ wf_left = wf(0)
+ J_m_left = sum(J_m_bound(p) for p in wf_left)
+ J_z_left = sum(J_z_bound(p) for p in wf_left)
+#HIDDEN_END_bind
+
+ plot_currents(syst, [
+ (r'$J_{\mathbf{m}_i}$ (from left)', J_m_left),
+ (r'$J_{σ_z}$ (from left)', J_z_left),
+ ])
+
+
+if __name__ == '__main__':
+ main()
diff --git a/doc/source/tutorial/plot_graphene.py b/doc/source/tutorial/plot_graphene.py
index d9181734ba869a624d03617c565f0852c3c50792..e81e331b2bedf7de7241b857a50611a562d51316 100644
--- a/doc/source/tutorial/plot_graphene.py
+++ b/doc/source/tutorial/plot_graphene.py
@@ -1,4 +1,4 @@
-# Tutorial 2.7.1. 2D example: graphene quantum dot
+# Tutorial 2.8.1. 2D example: graphene quantum dot
# ================================================
#
# Physics background
diff --git a/doc/source/tutorial/plot_zincblende.py b/doc/source/tutorial/plot_zincblende.py
index 467f5767e2d5d2e0e4d08d168c6ba7b2de5f69f0..6b75a54e0df3d4b977ab1ec51a00960ef4f6f262 100644
--- a/doc/source/tutorial/plot_zincblende.py
+++ b/doc/source/tutorial/plot_zincblende.py
@@ -1,4 +1,4 @@
-# Tutorial 2.7.2. 3D example: zincblende structure
+# Tutorial 2.8.2. 3D example: zincblende structure
# ================================================
#
# Physical background
diff --git a/doc/source/tutorial/tutorial6.rst b/doc/source/tutorial/tutorial6.rst
index dbf267042eb65f06d614fe3a04d18594b5cd3ec2..c2c314acd3f9a89464ef2ea67d503fbe4d111b4e 100644
--- a/doc/source/tutorial/tutorial6.rst
+++ b/doc/source/tutorial/tutorial6.rst
@@ -1,286 +1,367 @@
-Plotting Kwant systems and data in various styles
--------------------------------------------------
-
-The plotting functionality of Kwant has been used extensively (through
-`~kwant.plotter.plot` and `~kwant.plotter.map`) in the previous tutorials. In
-addition to this basic use, `~kwant.plotter.plot` offers many options to change
-the plotting style extensively. It is the goal of this tutorial to show how
-these options can be used to achieve various very different objectives.
-
-2D example: graphene quantum dot
-................................
+Computing local quantities: densities and currents
+==================================================
+In the previous tutorials we have mainly concentrated on calculating *global*
+properties such as conductance and band structures. Often, however, insight can
+be gained from calculating *locally-defined* quantities, that is, quantities
+defined over individual sites or hoppings in your system. In the
+:ref:`closed-systems` tutorial we saw how we could visualize the density
+associated with the eigenstates of a system using `kwant.plotter.map`.
+
+In this tutorial we will see how we can calculate more general quantities than
+simple densities by studying spin transport in a system with a magnetic
+texture.
.. seealso::
The complete source code of this example can be found in
- :download:`tutorial/plot_graphene.py <../../../tutorial/plot_graphene.py>`
-
-We begin by first considering a circular graphene quantum dot (similar to what
-has been used in parts of the tutorial :ref:`tutorial-graphene`.) In contrast
-to previous examples, we will also use hoppings beyond next-nearest neighbors:
+ :download:`tutorial/magnetic_texture.py <../../../tutorial/magnetic_texture.py>`
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_makesyst
- :end-before: #HIDDEN_END_makesyst
-Note that adding hoppings hoppings to the `n`-th nearest neighbors can be
-simply done by passing `n` as an argument to
-`~kwant.lattice.Polyatomic.neighbors`. Also note that we use the method
-`~kwant.builder.Builder.eradicate_dangling` to get rid of single atoms sticking
-out of the shape. It is necessary to do so *before* adding the
-next-nearest-neighbor hopping [#]_.
+Introduction
+------------
+Our starting point will be the following spinful tight-binding model on
+a square lattice:
-Of course, the system can be plotted simply with default settings:
+.. math::
+ H = - \sum_{⟨ij⟩}\sum_{α} |iα⟩⟨jα|
+ + J \sum_{i}\sum_{αβ} \mathbf{m}_i⋅ \mathbf{σ}_{αβ} |iα⟩⟨iβ|,
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotsyst1
- :end-before: #HIDDEN_END_plotsyst1
+where latin indices run over sites, and greek indices run over spin. We can
+identify the first term as a nearest-neighbor hopping between like-spins, and
+the second as a term that couples spins on the same site. The second term acts
+like a magnetic field of strength :math:`J` that varies from site to site and
+that, on site :math:`i`, points in the direction of the unit vector
+:math:`\mathbf{m}_i`. :math:`\mathbf{σ}_{αβ}` is a vector of Pauli matrices.
+We shall take the following form for :math:`\mathbf{m}_i`:
-However, due to the richer structure of the lattice, this results in a rather
-busy plot:
+.. math::
+ \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}{δ},
-.. image:: ../images/plot_graphene_syst1.*
+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}`.
-A much clearer plot can be obtained by using different colors for both
-sublattices, and by having different line widths for different hoppings. This
-can be achieved by passing a function to the arguments of
-`~kwant.plotter.plot`, instead of a constant. For properties of sites, this
-must be a function taking one site as argument, for hoppings a function taking
-the start end end site of hopping as arguments:
+To define this model in Kwant we start as usual by defining functions that
+depend on the model parameters:
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotsyst2
- :end-before: #HIDDEN_END_plotsyst2
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_model
+ :end-before: #HIDDEN_END_model
-Note that since we are using an unfinalized Builder, a ``site`` is really an
-instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
-that carries the same information, but is much easier to interpret:
+and define our system as a square shape on a square lattice with two orbitals
+per site, with leads attached on the left and right:
-.. image:: ../images/plot_graphene_syst2.*
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_syst
+ :end-before: #HIDDEN_END_syst
-Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
-to plot *data* living on the system.
+Below is a plot of a projection of :math:`\mathbf{m}_i` onto the x-y plane
+inside the scattering region. The z component is shown by the color scale:
-As an example, we now compute the eigenstates of the graphene quantum dot and
-intend to plot the wave function probability in the quantum dot. For aesthetic
-reasons (the wave functions look a bit nicer), we restrict ourselves to
-nearest-neighbor hopping. Computing the wave functions is done in the usual
-way (note that for a large-scale system, one would probably want to use sparse
-linear algebra):
+.. image:: ../images/mag_field_direction.*
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata1
- :end-before: #HIDDEN_END_plotdata1
+We will now be interested in analyzing the form of the scattering states
+that originate from the left lead:
-In most cases, to plot the wave function probability, one wouldn't use
-`~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
-`n`-th wave function using it:
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_wavefunction
+ :end-before: #HIDDEN_END_wavefunction
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata2
- :end-before: #HIDDEN_END_plotdata2
+Local densities
+---------------
+If we were simulating a spinless system with only a single degree of freedom,
+then calculating the density on each site would be as simple as calculating the
+absolute square of the wavefunction like::
-This results in a standard pseudocolor plot, showing in this case (``n=225``) a
-graphene edge state, i.e. a wave function mostly localized at the zigzag edges
-of the quantum dot.
+ density = np.abs(psi)**2
-.. image:: ../images/plot_graphene_data1.*
+When there are multiple degrees of freedom per site, however, one has to be
+more careful. In the present case with two (spin) degrees of freedom per site
+one could calculate the per-site density like:
-However although in general preferable, `~kwant.plotter.map` has a few
-deficiencies for this small system: For example, there are a few distortions at
-the edge of the dot. (This cannot be avoided in the type of interpolation used
-in `~kwant.plotter.map`). However, we can also use `~kwant.plotter.plot` to
-achieve a similar, but smoother result.
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_ldos
+ :end-before: #HIDDEN_END_ldos
-For this note that `~kwant.plotter.plot` can also take an array of floats (or
-function returning floats) as value for the ``site_color`` argument (the same
-holds for the hoppings). Via the colormap specified in ``cmap`` these are mapped
-to color, just as `~kwant.plotter.map` does! In addition, we can also change
-the symbol shape depending on the sublattice. With a triangle pointing up and
-down on the respective sublattice, the symbols used by plot fill the space
-completely:
+With more than one degree of freedom per site we have more freedom as to what
+local quantities we can meaningfully compute. For example, we may wish to
+calculate the local z-projected spin density. We could calculate
+this in the following way:
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata3
- :end-before: #HIDDEN_END_plotdata3
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_lsdz
+ :end-before: #HIDDEN_END_lsdz
-Note that with ``hop_lw=0`` we deactivate plotting the hoppings (that would not
-serve any purpose here). Moreover, ``site_size=0.5`` guarantees that the two
-different types of triangles touch precisely: By default, `~kwant.plotter.plot`
-takes all sizes in units of the nearest-neighbor spacing. ``site_size=0.5``
-thus means half the distance between neighboring sites (and for the triangles
-this is interpreted as the radius of the inner circle).
+If we wanted instead to calculate the local y-projected spin density, we would
+need to use an even more complicated expression:
-Finally, note that since we are dealing with a finalized system now, a site `i`
-is represented by an integer. In order to obtain the original
-`~kwant.builder.Site`, ``syst.sites[i]`` can be used.
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_lsdy
+ :end-before: #HIDDEN_END_lsdy
-With this we arrive at
+The `kwant.operator` module aims to alleviate somewhat this tedious
+book-keeping by providing a simple interface for defining operators that act on
+wavefunctions. To calculate the above quantities we would use the
+`~kwant.operator.Density` operator like so:
-.. image:: ../images/plot_graphene_data2.*
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_lden
+ :end-before: #HIDDEN_END_lden
-with the same information as `~kwant.plotter.map`, but with a cleaner look.
+`~kwant.operator.Density` takes a `~kwant.system.System` as its first parameter
+as well as (optionally) a square matrix that defines the quantity that you wish
+to calculate per site. When an instance of a `~kwant.operator.Density` is then
+evaluated with a wavefunction, the quantity
-The way how data is presented of course influences what features of the data
-are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
-beyond pseudocolor-like plots. For example, we can represent the wave function
-probability using the symbols itself:
+.. math:: ρ_i = \mathbf{ψ}^†_i \mathbf{M} \mathbf{ψ}_i
-.. literalinclude:: plot_graphene.py
- :start-after: #HIDDEN_BEGIN_plotdata4
- :end-before: #HIDDEN_END_plotdata4
+is calculated for each site :math:`i`, where :math:`\mathbf{ψ}_{i}` is a vector
+consisting of the wavefunction components on that site and :math:`\mathbf{M}`
+is the square matrix referred to previously.
-Here, we choose the symbol size proportional to the wave function probability,
-while the site color is transparent to also allow for overlapping symbols to be
-visible. The hoppings are also plotted in order to show the underlying lattice.
+Below we can see colorplots of the above-calculated quantities. The array that
+is returned by evaluating a `~kwant.operator.Density` can be used directly with
+`kwant.plotter.map`:
-With this, we arrive at
+.. image:: ../images/spin_densities.*
-.. image:: ../images/plot_graphene_data3.*
-which shows the edge state nature of the wave function most clearly.
+.. specialnote:: Technical Details
-.. rubric:: Footnotes
+ Although we refer loosely to "densities" and "operators" above, a
+ `~kwant.operator.Density` actually represents a *collection* of linear
+ operators. This can be made clear by rewriting the above definition
+ of :math:`ρ_i` in the following way:
-.. [#] A dangling site is defined as having only one hopping connecting it to
- the rest. With next-nearest-neighbor hopping also all sites that are
- dangling with only nearest-neighbor hopping have more than one hopping.
+ .. math::
+ ρ_i = \sum_{αβ} ψ^*_{α} \mathcal{M}_{iαβ} ψ_{β}
-3D example: zincblende structure
-................................
+ where greek indices run over the degrees of freedom in the Hilbert space of
+ the scattering region and latin indices run over sites. We can this
+ identify :math:`\mathcal{M}_{iαβ}` as the components of a rank-3 tensor and can
+ represent them as a "vector of matrices":
-.. seealso::
- The complete source code of this example can be found in
- :download:`tutorial/plot_zincblende.py <../../../tutorial/plot_zincblende.py>`
+ .. math::
+ \mathcal{M} = \left[
+ \left(\begin{matrix}
+ \mathbf{M} & 0 & … \\
+ 0 & 0 & … \\
+ ⋮ & ⋮ & ⋱
+ \end{matrix}\right)
+ ,\
+ \left(\begin{matrix}
+ 0 & 0 & … \\
+ 0 & \mathbf{M} & … \\
+ ⋮ & ⋮ & ⋱
+ \end{matrix}\right)
+ , … \right]
+
+ where :math:`\mathbf{M}` is defined as in the main text, and the :math:`0`
+ are zero matrices of the same shape as :math:`\mathbf{M}`.
-Zincblende is a very common crystal structure of semiconductors. It is a
-face-centered cubic crystal with two inequivalent atoms in the unit cell
-(i.e. two different types of atoms, unlike diamond which has the same crystal
-structure, but two equivalent atoms per unit cell).
-It is very easily generated in Kwant with `kwant.lattice.general`:
+Local currents
+--------------
+`kwant.operator` also has a class `~kwant.operator.Current` for calculating
+local currents, analogously to the local "densities" described above. If
+one has defined a density via a matrix :math:`\mathbf{M}` and the above
+equation, then one can define a local current flowing from site :math:`b`
+to site :math:`a`:
+
+.. math:: J_{ab} = i \left(
+ \mathbf{ψ}^†_b (\mathbf{H}_{ab})^† \mathbf{M} \mathbf{ψ}_a
+ - \mathbf{ψ}^†_a \mathbf{M} \mathbf{H}_{ab} \mathbf{ψ}_b
+ \right),
-.. literalinclude:: plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_zincblende1
- :end-before: #HIDDEN_END_zincblende1
+where :math:`\mathbf{H}_{ab}` is the hopping matrix from site :math:`b` to site
+:math:`a`. For example, to calculate the local current and
+spin current:
+
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_current
+ :end-before: #HIDDEN_END_current
-Note how we keep references to the two different sublattices for later use.
+Evaluating a `~kwant.operator.Current` operator on a wavefunction returns a
+1D array of values that can be directly used with `kwant.plotter.current`:
-A three-dimensional structure is created as easily as in two dimensions, by
-using the `~kwant.lattice.PolyatomicLattice.shape`-functionality:
+.. image:: ../images/spin_currents.*
-.. literalinclude:: plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_zincblende2
- :end-before: #HIDDEN_END_zincblende2
+.. note::
-We restrict ourselves here to a simple cuboid, and do not bother to add real
-values for onsite and hopping energies, but only the placeholder ``None`` (in a
-real calculation, several atomic orbitals would have to be considered).
+ Evaluating a `~kwant.operator.Current` operator on a wavefunction
+ returns a 1D array of the same length as the number of hoppings in the
+ system, ordered in the same way as the edges in the system's graph.
-`~kwant.plotter.plot` can plot 3D systems just as easily as its two-dimensional
-counterparts:
+.. specialnote:: Technical Details
-.. literalinclude:: plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_plot1
- :end-before: #HIDDEN_END_plot1
+ Similarly to how we saw in the previous section that `~kwant.operator.Density`
+ can be thought of as a collection of operators, `~kwant.operator.Current`
+ can be defined in a similar way. Starting from the definition of a "density":
-resulting in
+ .. math:: ρ_a = \sum_{αβ} ψ^*_{α} \mathcal{M}_{aαβ} ψ_{β},
-.. image:: ../images/plot_zincblende_syst1.*
+ we can define *currents* :math:`J_{ab}` via the continuity equation:
-You might notice that the standard options for plotting are quite different in
-3D than in 2D. For example, by default hoppings are not printed, but sites are
-instead represented by little "balls" touching each other (which is achieved by
-a default ``site_size=0.5``). In fact, this style of plotting 3D shows quite
-decently the overall geometry of the system.
+ .. math:: \frac{∂ρ_a}{∂t} - \sum_{b} J_{ab} = 0
-When plotting into a window, the 3D plots can also be rotated and scaled
-arbitrarily, allowing for a good inspection of the geometry from all sides.
+ where the sum runs over sites :math:`b` neigboring site :math:`a`.
+ Plugging in the definition for :math:`ρ_a`, along with the Schrödinger
+ equation and the assumption that :math:`\mathcal{M}` is time independent,
+ gives:
-.. note::
+ .. math:: J_{ab} = \sum_{αβ}
+ ψ^*_α \left(i \sum_{γ}
+ \mathcal{H}^*_{abγα} \mathcal{M}_{aγβ}
+ - \mathcal{M}_{aαγ} \mathcal{H}_{abγβ}
+ \right) ψ_β,
+
+ where latin indices run over sites and greek indices run over the Hilbert
+ space degrees of freedom, and
+
+ .. math:: \mathcal{H}_{ab} = \left(\begin{matrix}
+ ⋱ & ⋮ & ⋮ & ⋮ & ⋰ \\
+ ⋯ & ⋱ & 0 & \mathbf{H}_{ab} & ⋯ \\
+ ⋯ & 0 & ⋱ & 0 & ⋯ \\
+ ⋯ & 0 & 0 & ⋱ & ⋯ \\
+ ⋰ & ⋮ & ⋮ & ⋮ & ⋱
+ \end{matrix}\right).
+
+ i.e. :math:`\mathcal{H}_{ab}` is a matrix that is zero everywhere
+ except on elements connecting *from* site :math:`b` *to* site :math:`a`,
+ where it is equal to the hopping matrix :math:`\mathbf{H}_{ab}` between
+ these two sites.
+
+ This allows us to identify the rank-4 quantity
+
+ .. math:: \mathcal{J}_{abαβ} = i \sum_{γ}
+ \mathcal{H}^*_{abγα} \mathcal{M}_{aγβ}
+ - \mathcal{M}_{aαγ} \mathcal{H}_{abγβ}
- Interactive 3D plots usually do not have the proper aspect ratio, but are a
- bit squashed. This is due to bugs in matplotlib's 3D plotting module that
- does not properly honor the corresponding arguments. By resizing the plot
- window however one can manually adjust the aspect ratio.
+ as the local current between connected sites.
-Also for 3D it is possible to customize the plot. For example, we
-can explicitly plot the hoppings as lines, and color sites differently
-depending on the sublattice:
+ The diagonal part of this quantity, :math:`\mathcal{J}_{aa}`,
+ represents the extent to which the density defined by :math:`\mathcal{M}_a`
+ is not conserved on site :math:`a`. It can be calculated using
+ `~kwant.operator.Source`, rather than `~kwant.operator.Current`, which
+ only computes the off-diagonal part.
-.. literalinclude:: plot_zincblende.py
- :start-after: #HIDDEN_BEGIN_plot2
- :end-before: #HIDDEN_END_plot2
-which results in a 3D plot that allows to interactively (when plotted
-in a window) explore the crystal structure:
+Spatially varying operators
+---------------------------
+The above examples are reasonably simple in the sense that the book-keeping
+required to manually calculate the various densities and currents is still
+manageable. Now we shall look at the case where we wish to calculate some
+projected spin currents, but where the spin projection axis varies from place
+to place. More specifically, we want to visualize the spin current along the
+direction of :math:`\mathbf{m}_i`, which changes continuously over the whole
+scattering region.
-.. image:: ../images/plot_zincblende_syst2.*
+Doing this is as simple as passing a *function* when instantiating
+the `~kwant.operator.Current`, instead of a constant matrix:
-Hence, a few lines of code using Kwant allow to explore all the different
-crystal lattices out there!
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_following
+ :end-before: #HIDDEN_END_following
+
+The function must take a `~kwant.builder.Site` as its first parameter,
+and may optionally take other parameters (i.e. it must have the same
+signature as a Hamiltonian onsite function), and must return the square
+matrix that defines the operator we wish to calculate.
.. note::
- - The 3D plots are in fact only *fake* 3D. For example, sites will always
- be plotted above hoppings (this is due to the limitations of matplotlib's
- 3d module)
- - Plotting hoppings in 3D is inherently much slower than plotting sites.
- Hence, this is not done by default.
+ In the above example we had to pass the extra parameters needed by the
+ ``following_operator`` function via the ``param`` keyword argument. In
+ general you must pass all the parameters needed by the Hamiltonian via
+ ``params`` (as you would when calling `~kwant.solvers.default.smatrix` or
+ `~kwant.solvers.default.wave_function`). In the previous examples,
+ however, we used the fact that the system hoppings do not depend on any
+ parameters (these are the only Hamiltonian elements required to calculate
+ currents) to avoid passing the system parameters for the sake of brevity.
+Using this we can see that the spin current is essentially oriented along
+the direction of :math:`m_i` in the present regime where the onsite term
+in the Hamiltonian is dominant:
-Interpolated density and current: QPC with disorder
-...................................................
+.. image:: ../images/spin_current_comparison.*
-.. seealso::
- The complete source code of this example can be found in
- :download:`tutorial/plot_qpc.py <../../../tutorial/plot_qpc.py>`
+.. note:: Although this example used exclusively `~kwant.operator.Current`,
+ you can do the same with `~kwant.operator.Density`.
-In the above examples we saw some useful methods for plotting systems where
-single-site resolution is required. Sometimes, however, having single-site
-precision is a hinderance, rather than a help, and looking at *averaged*
-quantities is more useful. This is particularly important in systems with a
-large number of sites, and systems that are discretizations of continuum
-models.
-Here we will show how to plot interpolated quantities using `kwant.plotter.map`
-and `kwant.plotter.current` using the example of a quantum point contact (QPC)
-with a perpendicular magnetic field and disorder:
+Defining operators over parts of a system
+-----------------------------------------
-.. literalinclude:: plot_qpc.py
- :start-after: #HIDDEN_BEGIN_syst
- :end-before: #HIDDEN_END_syst
+Another useful feature of `kwant.operator` is the ability to calculate
+operators over selected parts of a system. For example, we may wish to
+calculate the total density of states in a certain part
+of the system, or the current flowing through a cut in the system.
+We can do this selection when creating the operator by using the
+keyword parameter ``where``.
-.. image:: ../images/plot_qpc_syst.*
+Density of states in a circle
+*****************************
-Now we will compute the density of particles and current due to states
-originating in the left lead with energy 0.15.
+To calculate the density of states inside a circle of radius
+20 we can simply do:
-.. literalinclude:: plot_qpc.py
- :start-after: #HIDDEN_BEGIN_wf
- :end-before: #HIDDEN_END_wf
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_density_cut
+ :end-before: #HIDDEN_END_density_cut
-We can then plot the density using `~kwant.plotter.map`:
+.. literalinclude:: ../images/circle_dos.txt
-.. literalinclude:: plot_qpc.py
- :start-after: #HIDDEN_BEGIN_density
- :end-before: #HIDDEN_END_density
+note that we also provide ``sum=True``, which means that evaluating the
+operator on a wavefunction will produce a single scalar. This is semantically
+equivalent to providing ``sum=False`` (the default) and running ``numpy.sum``
+on the output.
-.. image:: ../images/plot_qpc_density.*
+Current flowing through a cut
+*****************************
-We pass ``method='linear'`` to ``map``, which produces a smoother plot than the
-default style. We see that density is concentrated on the edges of the sample,
-as we expect due to Hall effect induced by the perpendicular magnetic field.
+Below we calculate the probability current and z-projected spin current near
+the interfaces with the left and right leads.
-Plotting the current in the system will enable us to make even more sense
-of what is going on:
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_current_cut
+ :end-before: #HIDDEN_END_current_cut
-.. literalinclude:: plot_qpc.py
- :start-after: #HIDDEN_BEGIN_current
- :end-before: #HIDDEN_END_current
+.. literalinclude:: ../images/current_cut.txt
+
+We see that the probability current is conserved across the scattering region,
+but the z-projected spin current is not due to the fact that the Hamiltonian
+does not commute with :math:`σ_z` everywhere in the scattering region.
+
+.. note:: ``where`` can also be provided as a sequence of `~kwant.builder.Site`
+ or a sequence of hoppings (i.e. pairs of `~kwant.builder.Site`),
+ rather than a function.
+
+
+Advanced Topics
+---------------
+
+Using ``bind`` for speed
+************************
+In most of the above examples we only used each operator *once* after creating
+it. Often one will want to evaluate an operator with many different
+wavefunctions, for example with all scattering wavefunctions at a certain
+energy, but with the *same set of parameters*. In such cases it is best to tell
+the operator to pre-compute the onsite matrices and any necessary Hamiltonian
+elements using the given set of parameters, so that this work is not duplicated
+every time the operator is evaluated.
+
+This can be achieved with `~kwant.operator.Current.bind`:
+
+.. warning:: Take care that you do not use an operator that was bound to a
+ particular set of parameters with wavefunctions calculated with a
+ *different* set of parameters. This will almost certainly give
+ incorrect results.
-.. image:: ../images/plot_qpc_current.*
+.. literalinclude:: magnetic_texture.py
+ :start-after: #HIDDEN_BEGIN_bind
+ :end-before: #HIDDEN_END_bind
-We can now clearly see that current enters the system from the lower-left edge
-of the system (this matches our intuition, as we calculated the current for
-scattering states originating in the left-hand lead), and is backscattered by
-the restriction of the QPC in the centre.
+.. image:: ../images/bound_current.*
diff --git a/doc/source/tutorial/tutorial7.rst b/doc/source/tutorial/tutorial7.rst
index 3ac75ba8936c8214ed2488246e045ef09d965263..dbf267042eb65f06d614fe3a04d18594b5cd3ec2 100644
--- a/doc/source/tutorial/tutorial7.rst
+++ b/doc/source/tutorial/tutorial7.rst
@@ -1,289 +1,286 @@
-##############################################################
-Calculating spectral density with the kernel polynomial method
-##############################################################
-
-We have already seen in the ":ref:`closed-systems`" tutorial that we can use
-Kwant simply to build Hamiltonians, which we can then directly diagonalize
-using routines from Scipy.
-
-This already allows us to treat systems with a few thousand sites without too
-many problems. For larger systems one is often not so interested in the exact
-eigenenergies and eigenstates, but more in the *density of states*.
-
-The kernel polynomial method (KPM), is an algorithm to obtain a polynomial
-expansion of the density of states. It can also be used to calculate the
-spectral density of arbitrary operators. Kwant has an implementation of the
-KPM method that is based on the algorithms presented in Ref. [1]_.
-
-Roughly speaking, KPM approximates the density of states (or any other spectral
-density) by expanding the action of the Hamiltonian (and operator of interest)
-on a (small) set of *random vectors* as a sum of Chebyshev polynomials up to
-some order, and then averaging. The accuracy of the method can be tuned by
-modifying the order of the Chebyshev expansion and the number of random
-vectors. See notes on accuracy_ below for details.
+Plotting Kwant systems and data in various styles
+-------------------------------------------------
+
+The plotting functionality of Kwant has been used extensively (through
+`~kwant.plotter.plot` and `~kwant.plotter.map`) in the previous tutorials. In
+addition to this basic use, `~kwant.plotter.plot` offers many options to change
+the plotting style extensively. It is the goal of this tutorial to show how
+these options can be used to achieve various very different objectives.
+
+2D example: graphene quantum dot
+................................
.. seealso::
The complete source code of this example can be found in
- :download:`tutorial/kernel_polynomial_method.py <../../../tutorial/kernel_polynomial_method.py>`
+ :download:`tutorial/plot_graphene.py <../../../tutorial/plot_graphene.py>`
+
+We begin by first considering a circular graphene quantum dot (similar to what
+has been used in parts of the tutorial :ref:`tutorial-graphene`.) In contrast
+to previous examples, we will also use hoppings beyond next-nearest neighbors:
-.. _accuracy:
-.. specialnote:: Performance and accuracy
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_makesyst
+ :end-before: #HIDDEN_END_makesyst
- The KPM method is especially well suited for large systems, and in the
- case when one is not interested in individual eigenvalues, but rather
- in obtaining an approximate spectral density.
+Note that adding hoppings hoppings to the `n`-th nearest neighbors can be
+simply done by passing `n` as an argument to
+`~kwant.lattice.Polyatomic.neighbors`. Also note that we use the method
+`~kwant.builder.Builder.eradicate_dangling` to get rid of single atoms sticking
+out of the shape. It is necessary to do so *before* adding the
+next-nearest-neighbor hopping [#]_.
- The accuracy in the energy resolution is dominated by the number of
- moments. The lowest accuracy is at the center of the spectrum, while
- slightly higher accuracy is obtained at the edges of the spectrum.
- If we use the KPM method (with the Jackson kernel, see Ref. [1]_) to
- describe a delta peak at the center of the spectrum, we will obtain a
- function similar to a Gaussian of width :math:`σ=πa/N`, where
- :math:`N` is the number of moments, and :math:`a` is the width of the
- spectrum.
+Of course, the system can be plotted simply with default settings:
- On the other hand, the random vectors will *explore* the range of the
- spectrum, and as the system gets bigger, the number of random vectors
- that are necessary to sample the whole spectrum reduces. Thus, a small
- number of random vectors is in general enough, and increasing this number
- will not result in a visible improvement of the approximation.
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotsyst1
+ :end-before: #HIDDEN_END_plotsyst1
+However, due to the richer structure of the lattice, this results in a rather
+busy plot:
-Introduction
-************
+.. image:: ../images/plot_graphene_syst1.*
-Our aim is to use the kernel polynomial method to obtain the spectral density
-:math:`ρ_A(E)`, as a function of the energy :math:`E`, of some Hilbert space
-operator :math:`A`. We define
+A much clearer plot can be obtained by using different colors for both
+sublattices, and by having different line widths for different hoppings. This
+can be achieved by passing a function to the arguments of
+`~kwant.plotter.plot`, instead of a constant. For properties of sites, this
+must be a function taking one site as argument, for hoppings a function taking
+the start end end site of hopping as arguments:
-.. math::
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotsyst2
+ :end-before: #HIDDEN_END_plotsyst2
- ρ_A(E) = ρ(E) A(E),
+Note that since we are using an unfinalized Builder, a ``site`` is really an
+instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
+that carries the same information, but is much easier to interpret:
-where :math:`A(E)` is the expectation value of :math:`A` for all the
-eigenstates of the Hamiltonian with energy :math:`E`, and the density of
-states is
+.. image:: ../images/plot_graphene_syst2.*
-.. math::
+Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
+to plot *data* living on the system.
- ρ(E) = \frac{1}{D} \sum_{k=0}^{D-1} δ(E-E_k),
+As an example, we now compute the eigenstates of the graphene quantum dot and
+intend to plot the wave function probability in the quantum dot. For aesthetic
+reasons (the wave functions look a bit nicer), we restrict ourselves to
+nearest-neighbor hopping. Computing the wave functions is done in the usual
+way (note that for a large-scale system, one would probably want to use sparse
+linear algebra):
-:math:`D` being the Hilbert space dimension, and :math:`E_k` the eigenvalues.
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotdata1
+ :end-before: #HIDDEN_END_plotdata1
-In the special case when :math:`A` is the identity, then :math:`ρ_A(E)` is
-simply :math:`ρ(E)`, the density of states.
+In most cases, to plot the wave function probability, one wouldn't use
+`~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
+`n`-th wave function using it:
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotdata2
+ :end-before: #HIDDEN_END_plotdata2
-Calculating the density of states
-*********************************
+This results in a standard pseudocolor plot, showing in this case (``n=225``) a
+graphene edge state, i.e. a wave function mostly localized at the zigzag edges
+of the quantum dot.
-In the following example, we will use the KPM implementation in Kwant
-to obtain the density of states of a graphene disk.
+.. image:: ../images/plot_graphene_data1.*
-We start by importing kwant and defining our system.
+However although in general preferable, `~kwant.plotter.map` has a few
+deficiencies for this small system: For example, there are a few distortions at
+the edge of the dot. (This cannot be avoided in the type of interpolation used
+in `~kwant.plotter.map`). However, we can also use `~kwant.plotter.plot` to
+achieve a similar, but smoother result.
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_sys1
- :end-before: #HIDDEN_END_sys1
+For this note that `~kwant.plotter.plot` can also take an array of floats (or
+function returning floats) as value for the ``site_color`` argument (the same
+holds for the hoppings). Via the colormap specified in ``cmap`` these are mapped
+to color, just as `~kwant.plotter.map` does! In addition, we can also change
+the symbol shape depending on the sublattice. With a triangle pointing up and
+down on the respective sublattice, the symbols used by plot fill the space
+completely:
-After making a system we can then create a `~kwant.kpm.SpectralDensity`
-object that represents the density of states for this system.
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotdata3
+ :end-before: #HIDDEN_END_plotdata3
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm1
- :end-before: #HIDDEN_END_kpm1
+Note that with ``hop_lw=0`` we deactivate plotting the hoppings (that would not
+serve any purpose here). Moreover, ``site_size=0.5`` guarantees that the two
+different types of triangles touch precisely: By default, `~kwant.plotter.plot`
+takes all sizes in units of the nearest-neighbor spacing. ``site_size=0.5``
+thus means half the distance between neighboring sites (and for the triangles
+this is interpreted as the radius of the inner circle).
-The `~kwant.kpm.SpectralDensity` can then be called like a function to obtain a
-sequence of energies in the spectrum of the Hamiltonian, and the corresponding
-density of states at these energies.
+Finally, note that since we are dealing with a finalized system now, a site `i`
+is represented by an integer. In order to obtain the original
+`~kwant.builder.Site`, ``syst.sites[i]`` can be used.
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm2
- :end-before: #HIDDEN_END_kpm2
+With this we arrive at
-When called with no arguments, an optimal set of energies is chosen (these are
-not evenly distributed over the spectrum, see Ref. [1]_ for details), however
-it is also possible to provide an explicit sequence of energies at which to
-evaluate the density of states.
+.. image:: ../images/plot_graphene_data2.*
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_kpm3
- :end-before: #HIDDEN_END_kpm3
+with the same information as `~kwant.plotter.map`, but with a cleaner look.
-.. image:: ../images/kpm_dos.*
+The way how data is presented of course influences what features of the data
+are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
+beyond pseudocolor-like plots. For example, we can represent the wave function
+probability using the symbols itself:
-In addition to being called like functions, `~kwant.kpm.SpectralDensity`
-objects also have a method `~kwant.kpm.SpectralDensity.average` which can be
-used to integrate the density of states against some distribution function over
-the whole spectrum. If no distribution function is specified, then the uniform
-distribution is used:
+.. literalinclude:: plot_graphene.py
+ :start-after: #HIDDEN_BEGIN_plotdata4
+ :end-before: #HIDDEN_END_plotdata4
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_av1
- :end-before: #HIDDEN_END_av1
+Here, we choose the symbol size proportional to the wave function probability,
+while the site color is transparent to also allow for overlapping symbols to be
+visible. The hoppings are also plotted in order to show the underlying lattice.
-.. literalinclude:: ../images/kpm_normalization.txt
+With this, we arrive at
-We see that the integral of the density of states is normalized to 1. If
-we wish to calculate, say, the average number of states populated in
-equilibrium, then we should integrate with respect to a Fermi-Dirac
-distribution and multiply by the total number of available states in
-the system:
+.. image:: ../images/plot_graphene_data3.*
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_av2
- :end-before: #HIDDEN_END_av2
+which shows the edge state nature of the wave function most clearly.
-.. literalinclude:: ../images/kpm_total_states.txt
+.. rubric:: Footnotes
-.. specialnote:: Stability and performance: spectral bounds
+.. [#] A dangling site is defined as having only one hopping connecting it to
+ the rest. With next-nearest-neighbor hopping also all sites that are
+ dangling with only nearest-neighbor hopping have more than one hopping.
- The KPM method internally rescales the spectrum of the Hamiltonian to the
- interval ``(-1, 1)`` (see Ref [1]_ for details), which requires calculating
- the boundaries of the spectrum (using ``scipy.sparse.linalg.eigsh``). This
- can be very costly for large systems, so it is possible to pass this
- explicitly as via the ``bounds`` parameter when instantiating the
- `~kwant.kpm.SpectralDensity` (see the class documentation for details).
+3D example: zincblende structure
+................................
- Additionally, `~kwant.kpm.SpectralDensity` accepts a parameter ``epsilon``,
- which ensures that the rescaled Hamiltonian (used internally), always has a
- spectrum strictly contained in the interval ``(-1, 1)``. If bounds are not
- provided then the tolerance on the bounds calculated with
- ``scipy.sparse.linalg.eigsh`` is set to ``epsilon/2``.
+.. seealso::
+ The complete source code of this example can be found in
+ :download:`tutorial/plot_zincblende.py <../../../tutorial/plot_zincblende.py>`
+Zincblende is a very common crystal structure of semiconductors. It is a
+face-centered cubic crystal with two inequivalent atoms in the unit cell
+(i.e. two different types of atoms, unlike diamond which has the same crystal
+structure, but two equivalent atoms per unit cell).
-Increasing the accuracy of the approximation
-********************************************
+It is very easily generated in Kwant with `kwant.lattice.general`:
-`~kwant.kpm.SpectralDensity` has two methods for increasing the accuracy
-of the method, each of which offers different levels of control over what
-exactly is changed.
+.. literalinclude:: plot_zincblende.py
+ :start-after: #HIDDEN_BEGIN_zincblende1
+ :end-before: #HIDDEN_END_zincblende1
-The simplest way to obtain a more accurate solution is to use the
-``add_moments`` method:
+Note how we keep references to the two different sublattices for later use.
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_acc1
- :end-before: #HIDDEN_END_acc1
+A three-dimensional structure is created as easily as in two dimensions, by
+using the `~kwant.lattice.PolyatomicLattice.shape`-functionality:
-This will update the number of calculated moments and also the
-number of sampling points such that the maximum distance between successive
-energy points is ``energy_resolution`` (see notes on accuracy_).
+.. literalinclude:: plot_zincblende.py
+ :start-after: #HIDDEN_BEGIN_zincblende2
+ :end-before: #HIDDEN_END_zincblende2
-.. image:: ../images/kpm_dos_acc.*
+We restrict ourselves here to a simple cuboid, and do not bother to add real
+values for onsite and hopping energies, but only the placeholder ``None`` (in a
+real calculation, several atomic orbitals would have to be considered).
-Alternatively, you can directly increase the number of moments
-with ``add_moments``, or the number of random vectors with ``add_vectors``.
+`~kwant.plotter.plot` can plot 3D systems just as easily as its two-dimensional
+counterparts:
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_acc2
- :end-before: #HIDDEN_END_acc2
+.. literalinclude:: plot_zincblende.py
+ :start-after: #HIDDEN_BEGIN_plot1
+ :end-before: #HIDDEN_END_plot1
-.. image:: ../images/kpm_dos_r.*
+resulting in
+.. image:: ../images/plot_zincblende_syst1.*
-.. _operator_spectral_density:
+You might notice that the standard options for plotting are quite different in
+3D than in 2D. For example, by default hoppings are not printed, but sites are
+instead represented by little "balls" touching each other (which is achieved by
+a default ``site_size=0.5``). In fact, this style of plotting 3D shows quite
+decently the overall geometry of the system.
-Calculating the spectral density of an operator
-***********************************************
+When plotting into a window, the 3D plots can also be rotated and scaled
+arbitrarily, allowing for a good inspection of the geometry from all sides.
-Above, we saw how to calculate the density of states by creating a
-`~kwant.kpm.SpectralDensity` and passing it a finalized Kwant system.
-When instantiating a `~kwant.kpm.SpectralDensity` we may optionally
-supply an operator in addition to the system. In this case it is
-the spectral density of the given operator that is calculated.
+.. note::
-`~kwant.kpm.SpectralDensity` accepts the operators in a few formats:
+ Interactive 3D plots usually do not have the proper aspect ratio, but are a
+ bit squashed. This is due to bugs in matplotlib's 3D plotting module that
+ does not properly honor the corresponding arguments. By resizing the plot
+ window however one can manually adjust the aspect ratio.
-* *explicit matrices* (numpy array of scipy sparse matrices will work)
-* *operators* from `kwant.operator`
+Also for 3D it is possible to customize the plot. For example, we
+can explicitly plot the hoppings as lines, and color sites differently
+depending on the sublattice:
-If an explicit matrix is provided then it must have the same
-shape as the system Hamiltonian.
+.. literalinclude:: plot_zincblende.py
+ :start-after: #HIDDEN_BEGIN_plot2
+ :end-before: #HIDDEN_END_plot2
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op1
- :end-before: #HIDDEN_END_op1
+which results in a 3D plot that allows to interactively (when plotted
+in a window) explore the crystal structure:
+.. image:: ../images/plot_zincblende_syst2.*
-Or, to do the same calculation using `kwant.operator.Density`:
+Hence, a few lines of code using Kwant allow to explore all the different
+crystal lattices out there!
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op2
- :end-before: #HIDDEN_END_op2
+.. note::
-Using operators from `kwant.operator` allows us to calculate quantities
-such as the *local* density of states by telling the operator not to
-sum over all the sites of the system:
+ - The 3D plots are in fact only *fake* 3D. For example, sites will always
+ be plotted above hoppings (this is due to the limitations of matplotlib's
+ 3d module)
+ - Plotting hoppings in 3D is inherently much slower than plotting sites.
+ Hence, this is not done by default.
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op3
- :end-before: #HIDDEN_END_op3
-`~kwant.kpm.SpectralDensity` will properly handle this vector output,
-which allows us to plot the local density of states at different
-point in the spectrum:
+Interpolated density and current: QPC with disorder
+...................................................
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_op4
- :end-before: #HIDDEN_END_op4
+.. seealso::
+ The complete source code of this example can be found in
+ :download:`tutorial/plot_qpc.py <../../../tutorial/plot_qpc.py>`
-.. image:: ../images/kpm_ldos.*
+In the above examples we saw some useful methods for plotting systems where
+single-site resolution is required. Sometimes, however, having single-site
+precision is a hinderance, rather than a help, and looking at *averaged*
+quantities is more useful. This is particularly important in systems with a
+large number of sites, and systems that are discretizations of continuum
+models.
-This nicely illustrates the edge states of the graphene dot at zero
-energy, and the bulk states at higher energy.
+Here we will show how to plot interpolated quantities using `kwant.plotter.map`
+and `kwant.plotter.current` using the example of a quantum point contact (QPC)
+with a perpendicular magnetic field and disorder:
+.. literalinclude:: plot_qpc.py
+ :start-after: #HIDDEN_BEGIN_syst
+ :end-before: #HIDDEN_END_syst
-Advanced topics
-***************
+.. image:: ../images/plot_qpc_syst.*
-Custom distributions for random vectors
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-By default `~kwant.kpm.SpectralDensity` will use random vectors
-whose components are unit complex numbers with phases drawn
-from a uniform distribution. There are several reasons why you may
-wish to make a different choice of distribution for your random vectors,
-for example to enforce certain symmetries or to only use real-valued vectors.
+Now we will compute the density of particles and current due to states
+originating in the left lead with energy 0.15.
-To change how the random vectors are generated, you need only specify a
-function that takes the dimension of the Hilbert space as a single parameter,
-and which returns a vector in that Hilbert space:
+.. literalinclude:: plot_qpc.py
+ :start-after: #HIDDEN_BEGIN_wf
+ :end-before: #HIDDEN_END_wf
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_fact1
- :end-before: #HIDDEN_END_fact1
+We can then plot the density using `~kwant.plotter.map`:
-Reproducible calculations
-^^^^^^^^^^^^^^^^^^^^^^^^^
-Because KPM internally uses random vectors, running the same calculation
-twice will not give bit-for-bit the same result. However, similarly to
-the funcions in `~kwant.rmt`, the random number generator can be directly
-manipulated by passing a value to the ``rng`` parameter of
-`~kwant.kpm.SpectralDensity`. ``rng`` can itself be a random number generator,
-or it may simply be a seed to pass to the numpy random number generator
-(that is used internally by default).
+.. literalinclude:: plot_qpc.py
+ :start-after: #HIDDEN_BEGIN_density
+ :end-before: #HIDDEN_END_density
-Defining operators as sesquilinear maps
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-`Above`__, we showed how `~kwant.kpm.SpectralDensity` can calculate the
-spectral density of operators, and how we can define operators by using
-`kwant.operator`. If you need even more flexibility,
-`~kwant.kpm.SpectralDensity` will also accept a *function* as its ``operator``
-parameter. This function must itself take two parameters, ``(bra, ket)`` and
-must return either a scalar or a one-dimensional array. In order to be
-meaningful the function must be a sesquilinear map, i.e. antilinear in its
-first argument, and linear in its second argument. Below, we compare two
-methods for computing the local density of states, one using
-`kwant.operator.Density`, and the other using a custom function.
+.. image:: ../images/plot_qpc_density.*
-.. literalinclude:: kernel_polynomial_method.py
- :start-after: #HIDDEN_BEGIN_blm
- :end-before: #HIDDEN_END_blm
+We pass ``method='linear'`` to ``map``, which produces a smoother plot than the
+default style. We see that density is concentrated on the edges of the sample,
+as we expect due to Hall effect induced by the perpendicular magnetic field.
-__ operator_spectral_density_
+Plotting the current in the system will enable us to make even more sense
+of what is going on:
+.. literalinclude:: plot_qpc.py
+ :start-after: #HIDDEN_BEGIN_current
+ :end-before: #HIDDEN_END_current
-.. rubric:: References
+.. image:: ../images/plot_qpc_current.*
-.. [1] `Rev. Mod. Phys., Vol. 78, No. 1 (2006)
- <https://arxiv.org/abs/cond-mat/0504627>`_.
+We can now clearly see that current enters the system from the lower-left edge
+of the system (this matches our intuition, as we calculated the current for
+scattering states originating in the left-hand lead), and is backscattered by
+the restriction of the QPC in the centre.
diff --git a/doc/source/tutorial/tutorial8.rst b/doc/source/tutorial/tutorial8.rst
index c210d9290fb158bb7f57f1a63bf853505b3debb1..f40d1920231ee7602e9d93c2f472a6e9cb023e22 100644
--- a/doc/source/tutorial/tutorial8.rst
+++ b/doc/source/tutorial/tutorial8.rst
@@ -1,323 +1,289 @@
-Discretizing continuous Hamiltonians
-------------------------------------
-
-Introduction
-............
-
-In ":ref:`tutorial_discretization_schrodinger`" we have learnt that Kwant works
-with tight-binding Hamiltonians. Often, however, one will start with a
-continuum model and will subsequently need to discretize it to arrive at a
-tight-binding model.
-Although discretizing a Hamiltonian is usually a simple
-process, it is tedious and repetitive. The situation is further exacerbated
-when one introduces additional on-site degrees of freedom, and tracking all
-the necessary terms becomes a chore.
-The `~kwant.continuum` sub-package aims to be a solution to this problem.
-It is a collection of tools for working with
-continuum models and for discretizing them into tight-binding models.
+##############################################################
+Calculating spectral density with the kernel polynomial method
+##############################################################
+
+We have already seen in the ":ref:`closed-systems`" tutorial that we can use
+Kwant simply to build Hamiltonians, which we can then directly diagonalize
+using routines from Scipy.
+
+This already allows us to treat systems with a few thousand sites without too
+many problems. For larger systems one is often not so interested in the exact
+eigenenergies and eigenstates, but more in the *density of states*.
+
+The kernel polynomial method (KPM), is an algorithm to obtain a polynomial
+expansion of the density of states. It can also be used to calculate the
+spectral density of arbitrary operators. Kwant has an implementation of the
+KPM method that is based on the algorithms presented in Ref. [1]_.
+
+Roughly speaking, KPM approximates the density of states (or any other spectral
+density) by expanding the action of the Hamiltonian (and operator of interest)
+on a (small) set of *random vectors* as a sum of Chebyshev polynomials up to
+some order, and then averaging. The accuracy of the method can be tuned by
+modifying the order of the Chebyshev expansion and the number of random
+vectors. See notes on accuracy_ below for details.
.. seealso::
- The complete source code of this tutorial can be found in
- :download:`tutorial/continuum_discretizer.py <../../../tutorial/continuum_discretizer.py>`
+ The complete source code of this example can be found in
+ :download:`tutorial/kernel_polynomial_method.py <../../../tutorial/kernel_polynomial_method.py>`
+.. _accuracy:
+.. specialnote:: Performance and accuracy
-.. _tutorial_discretizer_introduction:
+ The KPM method is especially well suited for large systems, and in the
+ case when one is not interested in individual eigenvalues, but rather
+ in obtaining an approximate spectral density.
-Discretizing by hand
-....................
+ The accuracy in the energy resolution is dominated by the number of
+ moments. The lowest accuracy is at the center of the spectrum, while
+ slightly higher accuracy is obtained at the edges of the spectrum.
+ If we use the KPM method (with the Jackson kernel, see Ref. [1]_) to
+ describe a delta peak at the center of the spectrum, we will obtain a
+ function similar to a Gaussian of width :math:`σ=πa/N`, where
+ :math:`N` is the number of moments, and :math:`a` is the width of the
+ spectrum.
-As an example, let us consider the following continuum Schrödinger equation
-for a semiconducting heterostructure (using the effective mass approximation):
+ On the other hand, the random vectors will *explore* the range of the
+ spectrum, and as the system gets bigger, the number of random vectors
+ that are necessary to sample the whole spectrum reduces. Thus, a small
+ number of random vectors is in general enough, and increasing this number
+ will not result in a visible improvement of the approximation.
-.. math::
- \left( k_x \frac{\hbar^2}{2 m(x)} k_x \right) \psi(x) = E \, \psi(x).
+Introduction
+************
-Replacing the momenta by their corresponding differential operators
+Our aim is to use the kernel polynomial method to obtain the spectral density
+:math:`ρ_A(E)`, as a function of the energy :math:`E`, of some Hilbert space
+operator :math:`A`. We define
.. math::
- k_\alpha = -i \partial_\alpha,
-
-for :math:`\alpha = x, y` or :math:`z`, and discretizing on a regular lattice of
-points with spacing :math:`a`, we obtain the tight-binding model
-.. math::
+ ρ_A(E) = ρ(E) A(E),
- H = - \frac{1}{a^2} \sum_i A\left(x+\frac{a}{2}\right)
- \big(\ket{i}\bra{i+1} + h.c.\big)
- + \frac{1}{a^2} \sum_i
- \left( A\left(x+\frac{a}{2}\right) + A\left(x-\frac{a}{2}\right)\right)
- \ket{i} \bra{i},
-
-with :math:`A(x) = \frac{\hbar^2}{2 m(x)}`.
-
-Using `~kwant.continuum.discretize` to obtain a template
-........................................................
-
-The function `kwant.continuum.discretize` takes a symbolic Hamiltonian and
-turns it into a `~kwant.builder.Builder` instance with appropriate spatial
-symmetry that serves as a template.
-(We will see how to use the template to build systems with a particular
-shape later).
-
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_symbolic_discretization
- :end-before: #HIDDEN_END_symbolic_discretization
-
-It is worth noting that ``discretize`` treats ``k_x`` and ``x`` as
-non-commuting operators, and so their order is preserved during the
-discretization process.
-
-Setting the ``verbose`` parameter to ``True`` prints extra information about the
-onsite and hopping functions assigned to the ``Builder`` produced
-by ``discretize``:
-
-.. literalinclude:: ../images/discretizer_intro_verbose.txt
-
-.. specialnote:: Technical details
-
- - ``kwant.continuum`` uses ``sympy`` internally to handle symbolic
- expressions. Strings are converted using `kwant.continuum.sympify`,
- which essentially applies some Kwant-specific rules (such as treating
- ``k_x`` and ``x`` as non-commutative) before calling ``sympy.sympify``
-
- - The builder returned by ``discretize`` will have an N-D
- translational symmetry, where ``N`` is the number of dimensions that were
- discretized. This is the case, even if there are expressions in the input
- (e.g. ``V(x, y)``) which in principle *may not* have this symmetry. When
- using the returned builder directly, or when using it as a template to
- construct systems with different/lower symmetry, it is important to
- ensure that any functional parameters passed to the system respect the
- symmetry of the system. Kwant provides no consistency check for this.
-
- - The discretization process consists of taking input
- :math:`H(k_x, k_y, k_z)`, multiplying it from the right by
- :math:`\psi(x, y, z)` and iteratively applying a second-order accurate
- central derivative approximation for every
- :math:`k_\alpha=-i\partial_\alpha`:
-
- .. math::
- \partial_\alpha \psi(\alpha) =
- \frac{1}{a} \left( \psi\left(\alpha + \frac{a}{2}\right)
- -\psi\left(\alpha - \frac{a}{2}\right)\right).
-
- This process is done separately for every summand in Hamiltonian.
- Once all symbols denoting operators are applied internal algorithm is
- calculating ``gcd`` for hoppings coming from each summand in order to
- find best possible approximation. Please see source code for details.
-
- - Instead of using ``discretize`` one can use
- `~kwant.continuum.discretize_symbolic` to obtain symbolic output.
- When working interactively in `Jupyter notebooks <https://jupyter.org/>`_
- it can be useful to use this to see a symbolic representation of
- the discretized Hamiltonian. This works best when combined with ``sympy``
- `Pretty Printing <http://docs.sympy.org/latest/tutorial/printing.html#setting-up-pretty-printing>`_.
-
- - The symbolic result of discretization obtained with
- ``discretize_symbolic`` can be converted into a
- builder using `~kwant.continuum.build_discretized`.
- This can be useful if one wants to alter the tight-binding Hamiltonian
- before building the system.
-
-
-Building a Kwant system from the template
-.........................................
-
-Let us now use the output of ``discretize`` as a template to
-build a system and plot some of its energy eigenstate. For this example the
-Hamiltonian will be
+where :math:`A(E)` is the expectation value of :math:`A` for all the
+eigenstates of the Hamiltonian with energy :math:`E`, and the density of
+states is
.. math::
- H = k_x^2 + k_y^2 + V(x, y),
-
-where :math:`V(x, y)` is some arbitrary potential.
-
-First, use ``discretize`` to obtain a
-builder that we will use as a template:
-
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_template
- :end-before: #HIDDEN_END_template
-
-We now use this system with the `~kwant.builder.Builder.fill`
-method of `~kwant.builder.Builder` to construct the system we
-want to investigate:
+ ρ(E) = \frac{1}{D} \sum_{k=0}^{D-1} δ(E-E_k),
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_fill
- :end-before: #HIDDEN_END_fill
+:math:`D` being the Hilbert space dimension, and :math:`E_k` the eigenvalues.
-After finalizing this system, we can plot one of the system's
-energy eigenstates:
+In the special case when :math:`A` is the identity, then :math:`ρ_A(E)` is
+simply :math:`ρ(E)`, the density of states.
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_plot_eigenstate
- :end-before: #HIDDEN_END_plot_eigenstate
-.. image:: ../images/discretizer_gs.*
+Calculating the density of states
+*********************************
-Note in the above that we provided the function ``V`` to
-``syst.hamiltonian_submatrix`` using ``params=dict(V=potential)``, rather than
-via ``args``.
+In the following example, we will use the KPM implementation in Kwant
+to obtain the density of states of a graphene disk.
-In addition, the function passed as ``V`` expects two input parameters ``x``
-and ``y``, the same as in the initial continuum Hamiltonian.
+We start by importing kwant and defining our system.
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_sys1
+ :end-before: #HIDDEN_END_sys1
-Models with more structure: Bernevig-Hughes-Zhang
-.................................................
+After making a system we can then create a `~kwant.kpm.SpectralDensity`
+object that represents the density of states for this system.
-When working with multi-band systems, like the Bernevig-Hughes-Zhang (BHZ)
-model [1]_ [2]_, one can provide matrix input to `~kwant.continuum.discretize`
-using ``identity`` and ``kron``. For example, the definition of the BHZ model can be
-written succinctly as:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_kpm1
+ :end-before: #HIDDEN_END_kpm1
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_define_qsh
- :end-before: #HIDDEN_END_define_qsh
+The `~kwant.kpm.SpectralDensity` can then be called like a function to obtain a
+sequence of energies in the spectrum of the Hamiltonian, and the corresponding
+density of states at these energies.
-We can then make a ribbon out of this template system:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_kpm2
+ :end-before: #HIDDEN_END_kpm2
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_define_qsh_build
- :end-before: #HIDDEN_END_define_qsh_build
+When called with no arguments, an optimal set of energies is chosen (these are
+not evenly distributed over the spectrum, see Ref. [1]_ for details), however
+it is also possible to provide an explicit sequence of energies at which to
+evaluate the density of states.
-and plot its dispersion using `kwant.plotter.bands`:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_kpm3
+ :end-before: #HIDDEN_END_kpm3
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_plot_qsh_band
- :end-before: #HIDDEN_END_plot_qsh_band
+.. image:: ../images/kpm_dos.*
-.. image:: ../images/discretizer_qsh_band.*
+In addition to being called like functions, `~kwant.kpm.SpectralDensity`
+objects also have a method `~kwant.kpm.SpectralDensity.average` which can be
+used to integrate the density of states against some distribution function over
+the whole spectrum. If no distribution function is specified, then the uniform
+distribution is used:
-In the above we see the edge states of the quantum spin Hall effect, which
-we can visualize using `kwant.plotter.map`:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_av1
+ :end-before: #HIDDEN_END_av1
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_plot_qsh_wf
- :end-before: #HIDDEN_END_plot_qsh_wf
+.. literalinclude:: ../images/kpm_normalization.txt
-.. image:: ../images/discretizer_qsh_wf.*
+We see that the integral of the density of states is normalized to 1. If
+we wish to calculate, say, the average number of states populated in
+equilibrium, then we should integrate with respect to a Fermi-Dirac
+distribution and multiply by the total number of available states in
+the system:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_av2
+ :end-before: #HIDDEN_END_av2
-Limitations of discretization
-.............................
+.. literalinclude:: ../images/kpm_total_states.txt
-It is important to remember that the discretization of a continuum
-model is an *approximation* that is only valid in the low-energy
-limit. For example, the quadratic continuum Hamiltonian
+.. specialnote:: Stability and performance: spectral bounds
-.. math::
-
- H_\textrm{continuous}(k_x) = \frac{\hbar^2}{2m}k_x^2
-
-
-and its discretized approximation
-
-.. math::
+ The KPM method internally rescales the spectrum of the Hamiltonian to the
+ interval ``(-1, 1)`` (see Ref [1]_ for details), which requires calculating
+ the boundaries of the spectrum (using ``scipy.sparse.linalg.eigsh``). This
+ can be very costly for large systems, so it is possible to pass this
+ explicitly as via the ``bounds`` parameter when instantiating the
+ `~kwant.kpm.SpectralDensity` (see the class documentation for details).
- H_\textrm{tight-binding}(k_x) = 2t \big(1 - \cos(k_x a)\big),
+ Additionally, `~kwant.kpm.SpectralDensity` accepts a parameter ``epsilon``,
+ which ensures that the rescaled Hamiltonian (used internally), always has a
+ spectrum strictly contained in the interval ``(-1, 1)``. If bounds are not
+ provided then the tolerance on the bounds calculated with
+ ``scipy.sparse.linalg.eigsh`` is set to ``epsilon/2``.
-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
-to how high in energy you need your tight-binding model to be valid.
+Increasing the accuracy of the approximation
+********************************************
-It is possible to set :math:`a` through the ``grid_spacing`` parameter
-to `~kwant.continuum.discretize`, as we will illustrate in the following
-example. Let us start from the continuum Hamiltonian
+`~kwant.kpm.SpectralDensity` has two methods for increasing the accuracy
+of the method, each of which offers different levels of control over what
+exactly is changed.
-.. math::
+The simplest way to obtain a more accurate solution is to use the
+``add_moments`` method:
- H(k) = k_x^2 \mathbb{1}_{2\times2} + α k_x \sigma_y.
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_acc1
+ :end-before: #HIDDEN_END_acc1
-We start by defining this model as a string and setting the value of the
-:math:`α` parameter:
+This will update the number of calculated moments and also the default
+number of sampling points such that the maximum distance between successive
+energy points is ``energy_resolution`` (see notes on accuracy_).
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_ls_def
- :end-before: #HIDDEN_END_ls_def
+.. image:: ../images/kpm_dos_acc.*
-Now we can use `kwant.continuum.lambdify` to obtain a function that computes
-:math:`H(k)`:
+Alternatively, you can directly increase the number of moments
+with ``add_moments``, or the number of random vectors with ``add_vectors``.
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_ls_hk_cont
- :end-before: #HIDDEN_END_ls_hk_cont
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_acc2
+ :end-before: #HIDDEN_END_acc2
-We can also construct a discretized approximation using
-`kwant.continuum.discretize`, in a similar manner to previous examples:
+.. image:: ../images/kpm_dos_r.*
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_ls_hk_tb
- :end-before: #HIDDEN_END_ls_hk_tb
-Below we can see the continuum and tight-binding dispersions for two
-different values of the discretization grid spacing :math:`a`:
+.. _operator_spectral_density:
-.. image:: ../images/discretizer_lattice_spacing.*
+Calculating the spectral density of an operator
+***********************************************
+Above, we saw how to calculate the density of states by creating a
+`~kwant.kpm.SpectralDensity` and passing it a finalized Kwant system.
+When instantiating a `~kwant.kpm.SpectralDensity` we may optionally
+supply an operator in addition to the system. In this case it is
+the spectral density of the given operator that is calculated.
-We clearly see that the smaller grid spacing is, the better we approximate
-the original continuous dispersion. It is also worth remembering that the
-Brillouin zone also scales with grid spacing: :math:`[-\frac{\pi}{a},
-\frac{\pi}{a}]`.
+`~kwant.kpm.SpectralDensity` accepts the operators in a few formats:
-.. specialnote:: Note
+* *explicit matrices* (numpy array of scipy sparse matrices will work)
+* *operators* from `kwant.operator`
- The use of `~kwant.wraparound.wraparound` makes it easy to extend this
- example to multidimensional band structures. ``wraparound`` replaces
- the translational symmetry in the ``x`` direction by a parameter ``k_x``.
+If an explicit matrix is provided then it must have the same
+shape as the system Hamiltonian.
-Advanced topics
-...............
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_op1
+ :end-before: #HIDDEN_END_op1
-The input to `kwant.continuum.discretize` and `kwant.continuum.lambdify` can be
-not only a ``string``, as we saw above, but also a ``sympy`` expression or
-a ``sympy`` matrix.
-This functionality will probably be mostly useful to people who
-are already experienced with ``sympy``.
+Or, to do the same calculation using `kwant.operator.Density`:
-It is possible to use ``identity`` (for identity matrix), ``kron`` (for Kronecker product), as well as Pauli matrices ``sigma_0``,
-``sigma_x``, ``sigma_y``, ``sigma_z`` in the input to
-`~kwant.continuum.lambdify` and `~kwant.continuum.discretize`, in order to simplify
-expressions involving matrices. Matrices can also be provided explicitly using
-square ``[]`` brackets. For example, all following expressions are equivalent:
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_op2
+ :end-before: #HIDDEN_END_op2
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_subs_1
- :end-before: #HIDDEN_END_subs_1
+Using operators from `kwant.operator` allows us to calculate quantities
+such as the *local* density of states by telling the operator not to
+sum over all the sites of the system:
-.. literalinclude:: ../images/discretizer_subs_1.txt
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_op3
+ :end-before: #HIDDEN_END_op3
-We can use the ``subs`` keyword parameter to substitute expressions
-and numerical values:
+`~kwant.kpm.SpectralDensity` will properly handle this vector output,
+which allows us to plot the local density of states at different
+point in the spectrum:
-.. literalinclude:: continuum_discretizer.py
- :start-after: #HIDDEN_BEGIN_subs_2
- :end-before: #HIDDEN_END_subs_2
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_op4
+ :end-before: #HIDDEN_END_op4
-.. literalinclude:: ../images/discretizer_subs_2.txt
+.. image:: ../images/kpm_ldos.*
-Symbolic expressions obtained in this way can be directly passed to all
-``discretizer`` functions.
+This nicely illustrates the edge states of the graphene dot at zero
+energy, and the bulk states at higher energy.
-.. specialnote:: Technical details
- Because of the way that ``sympy`` handles commutation relations all symbols
- representing position and momentum operators are set to be non commutative.
- This means that the order of momentum and position operators in the input
- expression is preserved. Note that it is not possible to define individual
- commutation relations within ``sympy``, even expressions such :math:`x k_y x`
- will not be simplified, even though mathematically :math:`[x, k_y] = 0`.
+Advanced topics
+***************
+
+Custom distributions for random vectors
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+By default `~kwant.kpm.SpectralDensity` will use random vectors
+whose components are unit complex numbers with phases drawn
+from a uniform distribution. There are several reasons why you may
+wish to make a different choice of distribution for your random vectors,
+for example to enforce certain symmetries or to only use real-valued vectors.
+
+To change how the random vectors are generated, you need only specify a
+function that takes the dimension of the Hilbert space as a single parameter,
+and which returns a vector in that Hilbert space:
+
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_fact1
+ :end-before: #HIDDEN_END_fact1
+
+Reproducible calculations
+^^^^^^^^^^^^^^^^^^^^^^^^^
+Because KPM internally uses random vectors, running the same calculation
+twice will not give bit-for-bit the same result. However, similarly to
+the funcions in `~kwant.rmt`, the random number generator can be directly
+manipulated by passing a value to the ``rng`` parameter of
+`~kwant.kpm.SpectralDensity`. ``rng`` can itself be a random number generator,
+or it may simply be a seed to pass to the numpy random number generator
+(that is used internally by default).
+
+Defining operators as sesquilinear maps
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+`Above`__, we showed how `~kwant.kpm.SpectralDensity` can calculate the
+spectral density of operators, and how we can define operators by using
+`kwant.operator`. If you need even more flexibility,
+`~kwant.kpm.SpectralDensity` will also accept a *function* as its ``operator``
+parameter. This function must itself take two parameters, ``(bra, ket)`` and
+must return either a scalar or a one-dimensional array. In order to be
+meaningful the function must be a sesquilinear map, i.e. antilinear in its
+first argument, and linear in its second argument. Below, we compare two
+methods for computing the local density of states, one using
+`kwant.operator.Density`, and the other using a custom function.
+
+.. literalinclude:: kernel_polynomial_method.py
+ :start-after: #HIDDEN_BEGIN_blm
+ :end-before: #HIDDEN_END_blm
+
+__ operator_spectral_density_
.. rubric:: References
-.. [1] `Science, 314, 1757 (2006)
- <https://arxiv.org/abs/cond-mat/0611399>`_.
-
-.. [2] `Phys. Rev. B 82, 045122 (2010)
- <https://arxiv.org/abs/1005.1682>`_.
+.. [1] `Rev. Mod. Phys., Vol. 78, No. 1 (2006)
+ <https://arxiv.org/abs/cond-mat/0504627>`_.
diff --git a/doc/source/tutorial/tutorial9.rst b/doc/source/tutorial/tutorial9.rst
new file mode 100644
index 0000000000000000000000000000000000000000..216f24b116f0d4036ac62b0ae77d37eccc302501
--- /dev/null
+++ b/doc/source/tutorial/tutorial9.rst
@@ -0,0 +1,318 @@
+Discretizing continuous Hamiltonians
+------------------------------------
+
+Introduction
+............
+
+In ":ref:`tutorial_discretization_schrodinger`" we have learnt that Kwant works
+with tight-binding Hamiltonians. Often, however, one will start with a
+continuum model and will subsequently need to discretize it to arrive at a
+tight-binding model.
+Although discretizing a Hamiltonian is usually a simple
+process, it is tedious and repetitive. The situation is further exacerbated
+when one introduces additional on-site degrees of freedom, and tracking all
+the necessary terms becomes a chore.
+The `~kwant.continuum` sub-package aims to be a solution to this problem.
+It is a collection of tools for working with
+continuum models and for discretizing them into tight-binding models.
+
+.. seealso::
+ The complete source code of this tutorial can be found in
+ :download:`tutorial/continuum_discretizer.py <../../../tutorial/continuum_discretizer.py>`
+
+
+.. _tutorial_discretizer_introduction:
+
+Discretizing by hand
+....................
+
+As an example, let us consider the following continuum Schrödinger equation
+for a semiconducting heterostructure (using the effective mass approximation):
+
+.. math::
+
+ \left( k_x \frac{\hbar^2}{2 m(x)} k_x \right) \psi(x) = E \, \psi(x).
+
+Replacing the momenta by their corresponding differential operators
+
+.. math::
+ k_\alpha = -i \partial_\alpha,
+
+for :math:`\alpha = x, y` or :math:`z`, and discretizing on a regular lattice of
+points with spacing :math:`a`, we obtain the tight-binding model
+
+.. math::
+
+ H = - \frac{1}{a^2} \sum_i A\left(x+\frac{a}{2}\right)
+ \big(\ket{i}\bra{i+1} + h.c.\big)
+ + \frac{1}{a^2} \sum_i
+ \left( A\left(x+\frac{a}{2}\right) + A\left(x-\frac{a}{2}\right)\right)
+ \ket{i} \bra{i},
+
+with :math:`A(x) = \frac{\hbar^2}{2 m(x)}`.
+
+Using `~kwant.continuum.discretize` to obtain a template
+........................................................
+
+The function `kwant.continuum.discretize` takes a symbolic Hamiltonian and
+turns it into a `~kwant.builder.Builder` instance with appropriate spatial
+symmetry that serves as a template.
+(We will see how to use the template to build systems with a particular
+shape later).
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_symbolic_discretization
+ :end-before: #HIDDEN_END_symbolic_discretization
+
+It is worth noting that ``discretize`` treats ``k_x`` and ``x`` as
+non-commuting operators, and so their order is preserved during the
+discretization process.
+
+Setting the ``verbose`` parameter to ``True`` prints extra information about the
+onsite and hopping functions assigned to the ``Builder`` produced
+by ``discretize``:
+
+.. literalinclude:: ../images/discretizer_intro_verbose.txt
+
+.. specialnote:: Technical details
+
+ - ``kwant.continuum`` uses ``sympy`` internally to handle symbolic
+ expressions. Strings are converted using `kwant.continuum.sympify`,
+ which essentially applies some Kwant-specific rules (such as treating
+ ``k_x`` and ``x`` as non-commutative) before calling ``sympy.sympify``
+
+ - The builder returned by ``discretize`` will have an N-D
+ translational symmetry, where ``N`` is the number of dimensions that were
+ discretized. This is the case, even if there are expressions in the input
+ (e.g. ``V(x, y)``) which in principle *may not* have this symmetry. When
+ using the returned builder directly, or when using it as a template to
+ construct systems with different/lower symmetry, it is important to
+ ensure that any functional parameters passed to the system respect the
+ symmetry of the system. Kwant provides no consistency check for this.
+
+ - The discretization process consists of taking input
+ :math:`H(k_x, k_y, k_z)`, multiplying it from the right by
+ :math:`\psi(x, y, z)` and iteratively applying a second-order accurate
+ central derivative approximation for every
+ :math:`k_\alpha=-i\partial_\alpha`:
+
+ .. math::
+ \partial_\alpha \psi(\alpha) =
+ \frac{1}{a} \left( \psi\left(\alpha + \frac{a}{2}\right)
+ -\psi\left(\alpha - \frac{a}{2}\right)\right).
+
+ This process is done separately for every summand in Hamiltonian.
+ Once all symbols denoting operators are applied internal algorithm is
+ calculating ``gcd`` for hoppings coming from each summand in order to
+ find best possible approximation. Please see source code for details.
+
+ - Instead of using ``discretize`` one can use
+ `~kwant.continuum.discretize_symbolic` to obtain symbolic output.
+ When working interactively in `Jupyter notebooks <https://jupyter.org/>`_
+ it can be useful to use this to see a symbolic representation of
+ the discretized Hamiltonian. This works best when combined with ``sympy``
+ `Pretty Printing <http://docs.sympy.org/latest/tutorial/printing.html#setting-up-pretty-printing>`_.
+
+ - The symbolic result of discretization obtained with
+ ``discretize_symbolic`` can be converted into a
+ builder using `~kwant.continuum.build_discretized`.
+ This can be useful if one wants to alter the tight-binding Hamiltonian
+ before building the system.
+
+
+Building a Kwant system from the template
+.........................................
+
+Let us now use the output of ``discretize`` as a template to
+build a system and plot some of its energy eigenstate. For this example the
+Hamiltonian will be
+
+.. math::
+
+ H = k_x^2 + k_y^2 + V(x, y),
+
+where :math:`V(x, y)` is some arbitrary potential.
+
+First, use ``discretize`` to obtain a
+builder that we will use as a template:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_template
+ :end-before: #HIDDEN_END_template
+
+We now use this system with the `~kwant.builder.Builder.fill`
+method of `~kwant.builder.Builder` to construct the system we
+want to investigate:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_fill
+ :end-before: #HIDDEN_END_fill
+
+After finalizing this system, we can plot one of the system's
+energy eigenstates:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_eigenstate
+ :end-before: #HIDDEN_END_plot_eigenstate
+
+.. image:: ../images/discretizer_gs.*
+
+Note in the above that we provided the function ``V`` to
+``syst.hamiltonian_submatrix`` using ``params=dict(V=potential)``, rather than
+via ``args``.
+
+In addition, the function passed as ``V`` expects two input parameters ``x``
+and ``y``, the same as in the initial continuum Hamiltonian.
+
+
+Models with more structure: Bernevig-Hughes-Zhang
+.................................................
+
+When working with multi-band systems, like the Bernevig-Hughes-Zhang (BHZ)
+model [1]_ [2]_, one can provide matrix input to `~kwant.continuum.discretize`
+using ``identity`` and ``kron``. For example, the definition of the BHZ model can be
+written succinctly as:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_define_qsh
+ :end-before: #HIDDEN_END_define_qsh
+
+We can then make a ribbon out of this template system:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_define_qsh_build
+ :end-before: #HIDDEN_END_define_qsh_build
+
+and plot its dispersion using `kwant.plotter.bands`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_qsh_band
+ :end-before: #HIDDEN_END_plot_qsh_band
+
+.. image:: ../images/discretizer_qsh_band.*
+
+In the above we see the edge states of the quantum spin Hall effect, which
+we can visualize using `kwant.plotter.map`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_plot_qsh_wf
+ :end-before: #HIDDEN_END_plot_qsh_wf
+
+.. image:: ../images/discretizer_qsh_wf.*
+
+
+Limitations of discretization
+.............................
+
+It is important to remember that the discretization of a continuum
+model is an *approximation* that is only valid in the low-energy
+limit. For example, the quadratic continuum Hamiltonian
+
+.. math::
+
+ H_\textrm{continuous}(k_x) = \frac{\hbar^2}{2m}k_x^2
+
+
+and its discretized approximation
+
+.. math::
+
+ H_\textrm{tight-binding}(k_x) = 2t \big(1 - \cos(k_x a)\big),
+
+
+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
+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
+to `~kwant.continuum.discretize`, as we will illustrate in the following
+example. Let us start from the continuum Hamiltonian
+
+.. math::
+
+ H(k) = k_x^2 \mathbb{1}_{2\times2} + α k_x \sigma_y.
+
+We start by defining this model as a string and setting the value of the
+:math:`α` parameter:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_def
+ :end-before: #HIDDEN_END_ls_def
+
+Now we can use `kwant.continuum.lambdify` to obtain a function that computes
+:math:`H(k)`:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_hk_cont
+ :end-before: #HIDDEN_END_ls_hk_cont
+
+We can also construct a discretized approximation using
+`kwant.continuum.discretize`, in a similar manner to previous examples:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_ls_hk_tb
+ :end-before: #HIDDEN_END_ls_hk_tb
+
+Below we can see the continuum and tight-binding dispersions for two
+different values of the discretization grid spacing :math:`a`:
+
+.. image:: ../images/discretizer_lattice_spacing.*
+
+
+We clearly see that the smaller grid spacing is, the better we approximate
+the original continuous dispersion. It is also worth remembering that the
+Brillouin zone also scales with grid spacing: :math:`[-\frac{\pi}{a},
+\frac{\pi}{a}]`.
+
+
+Advanced topics
+...............
+
+The input to `kwant.continuum.discretize` and `kwant.continuum.lambdify` can be
+not only a ``string``, as we saw above, but also a ``sympy`` expression or
+a ``sympy`` matrix.
+This functionality will probably be mostly useful to people who
+are already experienced with ``sympy``.
+
+
+It is possible to use ``identity`` (for identity matrix), ``kron`` (for Kronecker product), as well as Pauli matrices ``sigma_0``,
+``sigma_x``, ``sigma_y``, ``sigma_z`` in the input to
+`~kwant.continuum.lambdify` and `~kwant.continuum.discretize`, in order to simplify
+expressions involving matrices. Matrices can also be provided explicitly using
+square ``[]`` brackets. For example, all following expressions are equivalent:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_subs_1
+ :end-before: #HIDDEN_END_subs_1
+
+.. literalinclude:: ../images/discretizer_subs_1.txt
+
+We can use the ``substitutions`` keyword parameter to substitute expressions
+and numerical values:
+
+.. literalinclude:: continuum_discretizer.py
+ :start-after: #HIDDEN_BEGIN_subs_2
+ :end-before: #HIDDEN_END_subs_2
+
+.. literalinclude:: ../images/discretizer_subs_2.txt
+
+Symbolic expressions obtained in this way can be directly passed to all
+``discretizer`` functions.
+
+.. specialnote:: Technical details
+
+ Because of the way that ``sympy`` handles commutation relations all symbols
+ representing position and momentum operators are set to be non commutative.
+ This means that the order of momentum and position operators in the input
+ expression is preserved. Note that it is not possible to define individual
+ commutation relations within ``sympy``, even expressions such :math:`x k_y x`
+ will not be simplified, even though mathematically :math:`[x, k_y] = 0`.
+
+
+.. rubric:: References
+
+.. [1] `Science, 314, 1757 (2006)
+ <https://arxiv.org/abs/cond-mat/0611399>`_.
+
+.. [2] `Phys. Rev. B 82, 045122 (2010)
+ <https://arxiv.org/abs/1005.1682>`_.
diff --git a/kwant/operator.pyx b/kwant/operator.pyx
index 4ef1371eca0384ce0d9ab02a252d5862974db1ae..8fc66e81c5e856b274e2e422f7466e38fbe61952 100644
--- a/kwant/operator.pyx
+++ b/kwant/operator.pyx
@@ -815,7 +815,7 @@ cdef class Density(_LocalOperator):
cdef class Current(_LocalOperator):
- """An operator for calculating general currents.
+ r"""An operator for calculating general currents.
An instance of this class can be called like a function to evaluate the
expectation value with a wavefunction. See the documentation of the
@@ -850,9 +850,9 @@ cdef class Current(_LocalOperator):
represented by a square matrix :math:`M_i` associated with each site
:math:`i` and :math:`H_{ij}` is the hopping Hamiltonian from site :math:`j`
to site `i`, then an instance of this class represents the tensor
- :math:`J_{ijαβ}` which is equal to :math:`(H_{ij})^† M_i - M_i H_{ij}` when
- α and β are orbitals on sites :math:`i` and :math:`j` respectively, and
- zero otherwise.
+ :math:`J_{ijαβ}` which is equal to :math:`i\left[(H_{ij})^† M_i - M_i
+ H_{ij}\right]` when α and β are orbitals on sites :math:`i` and :math:`j`
+ respectively, and zero otherwise.
The tensor :math:`J_{ijαβ}` will also be referred to as :math:`Q_{nαβ}`,
where :math:`n` is the index of hopping :math:`(i, j)` in ``where``.