Commit 3659dad3 authored by Joseph Weston's avatar Joseph Weston

Merge branch 'll' into 'master'

Implement Landau level discretization.

See merge request kwant/kwant!295
parents 61d30d8f 4035315d
Pipeline #19044 passed with stages
in 8 minutes and 25 seconds
......@@ -12,3 +12,41 @@ the code are inserted directly into the document thanks to the magic of
`jupyter-sphinx <https://github.com/jupyter-widgets/jupyter-sphinx/>`_.
It has never been easier to get started contributing to Kwant by
helping us improve our documentation.
Discretization onto a Landau level basis
----------------------------------------
The Hamiltonian for a system infinite in at least two dimensions and with
a constant applied magnetic field may be expressed in a basis of Landau levels.
The momenta in the plane perpendicular to the magnetic field direction are
written in terms of the Landau level ladder operators:
.. math::
k_x = \sqrt{\frac{B}{2}} (a + a^\dagger) \quad\quad
k_y = i \sqrt{\frac{B}{2}} (a - a^\dagger)
The Hamiltonian is then expressed in terms of these ladder operators, which
allows for a straight-forward discretization in the basis of Landau levels,
provided that the basis is truncated after $N$ levels.
`kwant.continuum.discretize_landau` makes this procedure simple::
syst = kwant.continuum.discretize_landau("k_x**2 + k_y**2", N)
syst.finalized().hamiltonian_submatrix(params=dict(B=0.5))
`~kwant.continuum.discretize_landau` produces a regular Kwant Builder that
can be inspected or finalized as usual. 3D Hamiltonians for systems that
extend into the direction perpendicular to the magnetic field are also
possible::
template = kwant.continuum.discretize_landau("k_x**2 + k_y**2 + k_z**2 + V(z)", N)
This will produce a Builder with a single translational symmetry, which can be
directly finalized, or can be used as a template for e.g. a heterostructure stack
in the direction of the magnetic field::
def stack(site):
z, = site.pos
return 0 <= z < 10
template = kwant.continuum.discretize_landau("k_x**2 + k_y**2 + k_z**2 + V(z)", N)
syst = kwant.Builder()
syst.fill(template, stack, (0,))
......@@ -11,6 +11,7 @@ Discretizer
discretize
discretize_symbolic
build_discretized
discretize_landau
Symbolic helpers
----------------
......@@ -19,3 +20,11 @@ Symbolic helpers
sympify
lambdify
to_landau_basis
Other
-----
.. autosummary::
:toctree: generated/
LandauLattice
......@@ -217,6 +217,8 @@ In addition, the function passed as ``V`` expects two input parameters ``x``
and ``y``, the same as in the initial continuum Hamiltonian.
.. _discretize-bhz-model:
Models with more structure: Bernevig-Hughes-Zhang
.................................................
......
......@@ -12,4 +12,5 @@ Tutorial: learning Kwant through examples
plotting
kpm
discretize
magnetic_field
faq
Adding magnetic field
---------------------
Computing Landau levels in a harmonic oscillator basis
......................................................
When electrons move in an external magnetic field, their motion perpendicular
to the field direction is quantized into discrete Landau levels. Kwant implements
an efficient scheme for computing the Landau levels of arbitrary continuum
Hamiltonians. The general scheme revolves around rewriting the Hamiltonian in terms
of a basis of harmonic oscillator states [#]_, and is commonly illustrated in textbooks
for quadratic Hamiltonians.
To demonstrate the general scheme, let us consider a magnetic field oriented along
the :math:`z` direction :math:`\vec{B} = (0, 0, B)`, such that electron motion
in the :math:`xy` plane is Landau quantized. The magnetic field enters the Hamiltonian
through the kinetic momentum
.. math:: \hbar \vec{k} = - i \hbar \nabla + e\vec{A}(x, y).
In the symmetric gauge :math:`\vec{A}(x, y) = (B/2)[-y, x, 0]`, we introduce ladder
operators with the substitution
.. math::
k_x = \frac{1}{\sqrt{2} l_B} (a + a^\dagger), \quad \quad
k_y = \frac{1}{\sqrt{2} l_B} (a - a^\dagger),
with the magnetic length :math:`l_B = \sqrt{\hbar/eB}`. The ladder operators obey the
commutation relation
.. math:: [a, a^\dagger] = 1,
and define a quantum harmonic oscillator. We can thus write any electron continuum
Hamiltonian in terms of :math:`a` and :math:`a^\dagger`. Such a Hamiltonian has a
simple matrix representation in the eigenbasis of the number operator :math:`a^\dagger a`.
The eigenstates satisfy :math:`a^\dagger a | n \rangle = n | n \rangle` with the integer
Landau level index :math:`n \geq 0`, and in coordinate representation are proportional to
.. math::
\psi_n (x, y) = \left( \frac{\partial}{ \partial w} - \frac{w^*}{4 l_B^2} \right)
w^n e^{-|w|^2/4l_B^2},
with :math:`w = x + i y`. The matrix elements of the ladder operators are
.. math::
\langle n | a | m \rangle = \sqrt{m}~\delta_{n, m-1}, \quad \quad
\langle n | a^\dagger | m \rangle = \sqrt{m + 1}~\delta_{n, m+1}.
Truncating the basis to the first :math:`N` Landau levels allows us to approximate
the Hamiltonian as a discrete, finite matrix.
We can now formulate the algorithm that Kwant uses to compute Landau levels.
1. We take a generic continuum Hamiltonian, written in terms of the kinetic
momentum :math:`\vec{k}`. The Hamiltonian must be translationally
invariant along the directions perpendicular to the field direction.
2. We substitute the momenta perpendicular to the magnetic field with the ladder
operators :math:`a` and :math:`a^\dagger`.
3. We construct a `kwant.builder.Builder` using a special lattice which includes
the Landau level index :math:`n` as a degree of freedom on each site. The directions
normal to the field direction are not included in the builder, because they are
encoded in the Landau level index.
This procedure is automated with `kwant.continuum.discretize_landau`.
As an example, let us take the Bernevig-Hughes-Zhang model that we first considered in the
discretizer tutorial ":ref:`discretize-bhz-model`":
.. math::
C + M σ_0 \otimes σ_z + F(k_x^2 + k_y^2) σ_0 \otimes σ_z + D(k_x^2 + k_y^2) σ_0 \otimes σ_0
+ A k_x σ_z \otimes σ_x + A k_y σ_0 \otimes σ_y.
We can discretize this Hamiltonian in a basis of Landau levels as follows
.. jupyter-execute::
import numpy as np
import scipy.linalg
from matplotlib import pyplot
import kwant
import kwant.continuum
.. jupyter-execute:: boilerplate.py
:hide-code:
.. jupyter-execute::
hamiltonian = """
+ C * identity(4) + M * kron(sigma_0, sigma_z)
- F * (k_x**2 + k_y**2) * kron(sigma_0, sigma_z)
- D * (k_x**2 + k_y**2) * kron(sigma_0, sigma_0)
+ A * k_x * kron(sigma_z, sigma_x)
- A * k_y * kron(sigma_0, sigma_y)
"""
syst = kwant.continuum.discretize_landau(hamiltonian, N=10)
syst = syst.finalized()
We can then plot the spectrum of the system as a function of magnetic field, and
observe a crossing of Landau levels at finite magnetic field near zero energy,
characteristic of a quantum spin Hall insulator with band inversion.
.. jupyter-execute::
params = dict(A=3.645, F =-68.6, D=-51.2, M=-0.01, C=0)
b_values = np.linspace(0.0001, 0.0004, 200)
fig = kwant.plotter.spectrum(syst, ('B', b_values), params=params, show=False)
pyplot.ylim(-0.1, 0.2);
Comparing with tight-binding
============================
In the limit where fewer than one quantum of flux is threaded through a plaquette of
the discretization lattice we can compare the discretization in Landau levels with
a discretization in realspace.
.. jupyter-execute::
lat = kwant.lattice.square()
syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
def peierls(to_site, from_site, B):
y = from_site.tag[1]
return -1 * np.exp(-1j * B * y)
syst[(lat(0, j) for j in range(-19, 20))] = 4
syst[lat.neighbors()] = -1
syst[kwant.HoppingKind((1, 0), lat)] = peierls
syst = syst.finalized()
landau_syst = kwant.continuum.discretize_landau("k_x**2 + k_y**2", N=5)
landau_syst = landau_syst.finalized()
Here we plot the dispersion relation for the discretized ribbon and compare it
with the Landau levels shown as dashed lines.
.. jupyter-execute::
fig, ax = pyplot.subplots(1, 1)
ax.set_xlabel("momentum")
ax.set_ylabel("energy")
ax.set_ylim(0, 1)
params = dict(B=0.1)
kwant.plotter.bands(syst, ax=ax, params=params)
h = landau_syst.hamiltonian_submatrix(params=params)
for ev in scipy.linalg.eigvals(h):
ax.axhline(ev, linestyle='--')
The dispersion and the Landau levels diverge with increasing energy, because the real space
discretization of the ribbon gives a worse approximation to the dispersion at higher energies.
Discretizing 3D models
======================
Although the preceding examples have only included the plane perpendicular to the
magnetic field, the Landau level quantization also works if the direction
parallel to the field is included. In fact, although the system must be
translationally invariant in the plane perpendicular to the field, the system may
be arbitrary along the parallel direction. For example, it is therefore possible to
model a heterostructure and/or set up a scattering problem along the field direction.
Let's say that we wish to to model a heterostructure with a varying potential
:math:`V` along the direction of a magnetic field, :math:`z`, that includes
Zeeman splitting and Rashba spin-orbit coupling:
.. math::
\frac{\hbar^2}{2m}\sigma_0(k_x^2 + k_y^2 + k_z^2)
+ V(z)\sigma_0
+ \frac{\mu_B B}{2}\sigma_z
+ \hbar\alpha(\sigma_x k_y - \sigma_y k_x).
We can discretize this Hamiltonian in a basis of Landau levels as before:
.. jupyter-execute::
continuum_hamiltonian = """
(k_x**2 + k_y**2 + k_z**2) * sigma_0
+ V(z) * sigma_0
+ mu * B * sigma_z / 2
+ alpha * (sigma_x * k_y - sigma_y * k_x)
"""
template = kwant.continuum.discretize_landau(continuum_hamiltonian, N=10)
This creates a system with a single translational symmetry, along
the :math:`z` direction, which we can use as a template
to construct our heterostructure:
.. jupyter-execute::
def hetero_structure(site):
z, = site.pos
return 0 <= z < 10
def hetero_potential(z):
if z < 2:
return 0
elif z < 7:
return 0.5
else:
return 0.7
syst = kwant.Builder()
syst.fill(template, hetero_structure, (0,))
syst = syst.finalized()
params = dict(
B=0.5,
mu=0.2,
alpha=0.4,
V=hetero_potential,
)
syst.hamiltonian_submatrix(params=params);
.. rubric:: Footnotes
.. [#] `Wikipedia <https://en.wikipedia.org/wiki/Landau_quantization>`_ has
a nice introduction to Landau quantization.
\ No newline at end of file
......@@ -10,10 +10,12 @@ try:
from .discretizer import discretize, discretize_symbolic, build_discretized
from ._common import sympify, lambdify
from ._common import momentum_operators, position_operators
from .landau_levels import to_landau_basis, discretize_landau, LandauLattice
except ImportError as error:
msg = ("'kwant.continuum' is not available because one or more of its "
"dependencies is not installed.")
raise ImportError(msg) from error
__all__ = ['discretize', 'discretize_symbolic', 'build_discretized',
'to_landau_basis', 'discretize_landau', 'LandauLattice',
'sympify', 'lambdify', 'momentum_operators', 'position_operators']
This diff is collapsed.
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
# http://kwant-project.org/license. A list of Kwant authors can be found in
# the file AUTHORS.rst at the top-level directory of this distribution and at
# http://kwant-project.org/authors.
from math import sqrt
import numpy as np
import sympy
import pytest
import itertools
import kwant.builder
import kwant.lattice
from .._common import position_operators, momentum_operators, sympify
from ..landau_levels import (
to_landau_basis,
discretize_landau,
_ladder_term,
_evaluate_ladder_term,
ladder_lower,
ladder_raise,
LandauLattice,
)
x, y, z = position_operators
k_x, k_y, k_z = momentum_operators
B = sympy.symbols("B")
V = sympy.symbols("V", cls=sympy.Function)
a, a_dag = ladder_lower, ladder_raise
def test_ladder_term():
assert _ladder_term(a ** 2 * a_dag) == (-2, 1)
assert _ladder_term(a_dag ** 5 * a ** 3 * a_dag) == (5, -3, 1)
# non-ladder operators give a shift of 0
assert _ladder_term(B * a ** 2) == (0, -2)
def test_evaluate_ladder_term():
assert np.isclose(_evaluate_ladder_term((+1, -1), 1, +1), 1)
assert np.isclose(_evaluate_ladder_term((+1, -1), 2, +1), 2)
assert np.isclose(_evaluate_ladder_term((-2, +3, -2), 5, +1), 4 * 5 * 6 * sqrt(5))
# annihilating |0> is always 0
assert _evaluate_ladder_term((-1,), 0, +1) == 0
assert _evaluate_ladder_term((-2,), 1, +1) == 0
assert _evaluate_ladder_term((-3,), 1, +1) == 0
assert _evaluate_ladder_term((+3, -2), 1, +1) == 0
assert _evaluate_ladder_term((-3, -2, +3), 1, +1) == 0
# negative B swaps creation and annihilation operators
assert _evaluate_ladder_term((+1, -1), 2, +1) == _evaluate_ladder_term(
(-1, +1), 2, -1
)
assert _evaluate_ladder_term((-2, +3, -2), 5, +1) == _evaluate_ladder_term(
(+2, -3, +2), 5, -1
)
def test_to_landau_basis():
# test basic usage
ham, momenta, normal_coord = to_landau_basis("k_x**2 + k_y**2")
assert sympy.expand(ham) == abs(B) * a * a_dag + abs(B) * a_dag * a
assert momenta == (k_x, k_y)
assert normal_coord == z
# test that hamiltonian can be specified as a sympy expression
ham, momenta, normal_coord = to_landau_basis(sympify("k_x**2 + k_y**2"))
assert sympy.expand(ham) == abs(B) * a * a_dag + abs(B) * a_dag * a
assert momenta == (k_x, k_y)
assert normal_coord == z
# test that
ham, momenta, normal_coord = to_landau_basis("k_x**2 + k_y**2 + k_z**2 + V(z)")
assert sympy.expand(ham) == (
abs(B) * a * a_dag + abs(B) * a_dag * a + k_z ** 2 + V(z)
)
assert momenta == (k_x, k_y)
assert normal_coord == z
# test for momenta explicitly specified
ham, momenta, normal_coord = to_landau_basis(
"k_x**2 + k_y**2 + k_z**2 + k_x*k_y", momenta=("k_z", "k_y")
)
assert sympy.expand(ham) == (
abs(B) * a * a_dag
+ abs(B) * a_dag * a
+ k_x ** 2
+ sympy.I * sympy.sqrt(abs(B) / 2) * k_x * a
- sympy.I * sympy.sqrt(abs(B) / 2) * k_x * a_dag
)
assert momenta == (k_z, k_y)
assert normal_coord == x
def test_discretize_landau():
n_levels = 10
magnetic_field = 1 / 3 # a suitably arbitrary value
# Ensure that we can handle products of ladder operators by truncating
# several levels higher than the number of levels we actually want.
a = np.diag(np.sqrt(np.arange(1, n_levels + 5)), k=1)
a_dag = a.conjugate().transpose()
k_x = sqrt(magnetic_field / 2) * (a + a_dag)
k_y = 1j * sqrt(magnetic_field / 2) * (a - a_dag)
sigma_0 = np.eye(2)
sigma_x = np.array([[0, 1], [1, 0]])
sigma_y = np.array([[0, -1j], [1j, 0]])
# test that errors are raised on invalid input
with pytest.raises(ValueError):
discretize_landau("k_x**2 + k_y**2", N=0)
with pytest.raises(ValueError):
discretize_landau("k_x**2 + k_y**2", N=-1)
# test a basic Hamiltonian with no normal coordinate dependence
syst = discretize_landau("k_x**2 + k_y**2", N=n_levels)
lat = LandauLattice(1, norbs=1)
assert isinstance(syst.symmetry, kwant.builder.NoSymmetry)
syst = syst.finalized()
assert set(syst.sites) == {lat(0, j) for j in range(n_levels)}
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=0)), np.zeros((n_levels, n_levels))
)
should_be = k_x @ k_x + k_y @ k_y
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
should_be[:n_levels, :n_levels],
)
# test negative magnetic field
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=-magnetic_field)),
should_be[:n_levels, :n_levels],
)
# test hamiltonian with no onsite elements
syst = discretize_landau("k_x", N=n_levels)
syst = syst.finalized()
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
k_x[:n_levels, :n_levels],
)
# test a basic Hamiltonian with normal coordinate dependence
grid = 1 / 7 # a suitably arbitrary value
syst = discretize_landau(
"k_x**2 + k_y**2 + k_z**2 + V(z)", N=n_levels, grid_spacing=grid
)
assert isinstance(syst.symmetry, kwant.lattice.TranslationalSymmetry)
syst = syst.finalized()
zero_potential = syst.cell_hamiltonian(params=dict(B=magnetic_field, V=lambda z: 0))
with_potential = syst.cell_hamiltonian(params=dict(B=magnetic_field, V=lambda z: 1))
# extra +2/grid**2 from onsite kinetic energy
assert np.allclose(
zero_potential,
should_be[:n_levels, :n_levels] + (2 / grid ** 2) * np.eye(n_levels),
)
# V(z) just adds the same energy to each onsite
assert np.allclose(with_potential - zero_potential, np.eye(n_levels))
# hopping matrix does not exchange landau levels
assert np.allclose(
syst.inter_cell_hopping(params=dict(B=magnetic_field, V=lambda z: 0)),
-np.eye(n_levels) / grid ** 2,
)
# test a Hamiltonian with mixing between Landau levels
# and spatial degrees of freedom.
syst = discretize_landau("k_x**2 + k_y**2 + k_x*k_z", N=n_levels)
syst = syst.finalized()
assert np.allclose(
syst.inter_cell_hopping(params=dict(B=magnetic_field)),
-1j * k_x[:n_levels, :n_levels] / 2,
)
# test a Hamiltonian with extra degrees of freedom
syst = discretize_landau("sigma_0 * k_x**2 + sigma_x * k_y**2", N=n_levels)
syst = syst.finalized()
assert syst.sites[0].family.norbs == 2
should_be = np.kron(k_x @ k_x, sigma_0) + np.kron(k_y @ k_y, sigma_x)
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
should_be[: 2 * n_levels, : 2 * n_levels],
)
# test a linear Hamiltonian
syst = discretize_landau("sigma_y * k_x - sigma_x * k_y", N=n_levels)
syst = syst.finalized()
should_be = np.kron(k_x, sigma_y) - np.kron(k_y, sigma_x)
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
should_be[: 2 * n_levels, : 2 * n_levels],
)
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
syst.hamiltonian_submatrix(params=dict(B=-magnetic_field)),
)
def test_analytical_spectrum():
hamiltonian = """(k_x**2 + k_y**2) * sigma_0 +
alpha * (k_x * sigma_y - k_y * sigma_x) +
g * B * sigma_z"""
def exact_Es(n, B, alpha, g):
# See e.g. R. Winkler (2003), section 8.4.1
sign_B = np.sign(B)
B = np.abs(B)
Ep = 2*B*(n+1) - 0.5*np.sqrt((2*B - sign_B*2*g*B)**2 + 8*B*alpha**2*(n+1))
Em = 2*B*n + 0.5*np.sqrt((2*B - sign_B*2*g*B)**2 + 8*B*alpha**2*n)
return Ep, Em
N = 20
syst = discretize_landau(hamiltonian, N)
syst = syst.finalized()
params = dict(alpha=0.07, g=0.04)
for _ in range(5):
B = 0.01 + np.random.rand()*3
params['B'] = B
exact = [exact_Es(n, B, params['alpha'], params['g']) for n in range(N)]
# We only check the first N levels - the SOI couples adjacent levels,
# so the higher ones are not necessarily determined accurately in the
# discretization
exact = np.sort([energy for energies in exact for energy in energies])[:N]
ll_spect = np.linalg.eigvalsh(syst.hamiltonian_submatrix(params=params))[:len(exact)]
assert np.allclose(ll_spect, exact)
def test_fill():
def shape(site, lower, upper):
(z, )= site.pos
return lower <= z < upper
hamiltonian = "k_x**2 + k_y**2 + k_z**2"
N = 6
template = discretize_landau(hamiltonian, N)
syst = kwant.Builder()
width = 4
syst.fill(template, lambda site: shape(site, 0, width), (0, ));
correct_tags = [(coordinate, ll_index) for coordinate, ll_index
in itertools.product(range(width), range(N))]
syst_tags = [site.tag for site in syst.sites()]
assert len(syst_tags) == len(correct_tags)
assert all(tag in correct_tags for tag in syst_tags)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment