Tip

Run this Jupyter Notebook locally: Download Notebook

Tutorial Tiles

Tiles is a module to compute vectors that respect a certain symmetry, the symmetry of a superlattice. The only requirement for computations to make sense, is that the system is closed and has periodic boundaries, in which case it constitutes a superlattice.

These vectors are used to compute momentum-resolved spectral densities, or expectation values of any operator acting on the closed system.

The advantage of using tiles instad of local vectors, is that we can do Brillouin zone-unfolding; a procedure to compute the unfolded bands of superlattices.

Another benefit is that we can compute momentum-resolved spectral densities, even when the translation symmetry is broken or the system is disordered.

For example, if we would like to simulate the ARPES measurement of a non-pristine crystalline surface, the momentum components of a moire superlattice, then can use tiles to compute them.

from kpm_tools.tiles import (
    tile_random_kvectors, TileKVectors,
    momentum_to_ts
)
from kpm_tools.kpm_generator import SpectralDensityIterator

import matplotlib.pyplot, matplotlib.backends
import holoviews as hv
hv.extension()

import numpy as np

import kwant
from concurrent.futures import ProcessPoolExecutor
from kpm_tools.kpm_generator import _BaseKPM

We define globally the honeycomb lattice, in part to define a tile:

honeycomb = kwant.lattice.honeycomb(norbs=1, name=["a", "b"])

However, we don’t need to restrict the tiles to one unitcell of the lattice, it could span several unitcells, as long as it’s commensurate with the superlattice plus the periodic boundary conditions.

def make_superlattice(tags_1=[1, 0], tags_2=[0, 1], disorder=0):
    def disorder_onsite(site, disorder):
        return disorder * kwant.digest.uniform(bytes(str(site), "UTF8"))
    period_1 = honeycomb.prim_vecs.T @ tags_1
    period_2 = honeycomb.prim_vecs.T @ tags_2
    symm = kwant.TranslationalSymmetry(period_1, period_2)
    builder = kwant.Builder(symmetry=symm)

    builder[honeycomb.shape(lambda s: True, (0, 0))] = disorder_onsite

    builder[honeycomb.neighbors()] = 1

    return kwant.wraparound.wraparound(builder).finalized()

fsyst = make_superlattice(tags_1=[3, 0], tags_2=[0, 1], disorder=0.5)

kwant.plot(fsyst);

params = {'k_x': 0, 'k_y': 0, 'disorder': 0}

The bounds should be fixed, since we will iterate in reciprocal space

def init_kpm_params(fsyst, params):
    num_moments = 300
    kpm_params = {
        "num_moments": num_moments,
        "num_vectors": 8,
        "rng": 0,
    }

    bounds = _BaseKPM(
        fsyst,
        params=params,
        num_vectors=1,
        num_moments=2,
    ).bounds

    kpm_params["bounds"] = bounds


    energies = (
        kwant.kpm._chebyshev_nodes(2 * num_moments) * (bounds[1] - bounds[0]) / 2
        + (bounds[1] + bounds[0]) / 2
    )
    return kpm_params, energies

Define the k-path

def make_path_segment(p1, p2, num_k):
    return (
        p1[np.newaxis]
        + np.linspace(0, 1, num_k, endpoint=False)[:, np.newaxis]
        * (p2 - p1)[np.newaxis]
    )

def make_path(points, num_k):
    k_path = np.zeros((0, *points[0].shape))
    for idx in range(len(points)-1):
        p1 = points[idx]
        p2 = points[idx+1]
        k_path = np.concatenate([k_path, make_path_segment(p1, p2, num_k)], axis=0)
    return k_path

p_gamma = np.array([0 ,0])
p_m = np.array([1/2,-1/(2*np.sqrt(3))])
p_K = np.array([2/3, 0])
p_Kp = p_K + 2 * (p_m - p_K)

points = np.array([p_gamma, p_m, p_K, p_gamma, p_Kp]) * 2 * np.pi
num_k = 50

points = points

k_path = make_path(points, num_k)
num_k_full = len(k_path)

k_path_lengths = np.cumsum(np.append(0, np.linalg.norm(np.diff(k_path, axis=0), axis=1)))

xticks = [(k_path_lengths[num_k * 0], '$\Gamma$'),
          (k_path_lengths[num_k * 1], '$M$'),
          (k_path_lengths[num_k * 2], '$K$'),
          (k_path_lengths[num_k * 3], '$\Gamma$'),
          (k_path_lengths[num_k * 4-1], "$K'$"),
         ]

def density_in_k(k, tile):
    kvecs = tile_random_kvectors(k, tile, fsyst)

    # update boundary conditions to be periodic like `k`
    k_in_ts = momentum_to_ts(k, fsyst._wrapped_symmetry)
    params.update(dict(zip(fsyst._momentum_names, k_in_ts)))

    spectrum_k = SpectralDensityIterator(
        fsyst,
        params=params,
        vector_factory=kvecs,
        **kpm_params,
    )

    return spectrum_k.densities

Typical k-path plot

Here, we obtain the Brillouin zone-unfolded typical k-path for graphene or the honeycomb lattice.

Notice that the tile controls the translational symmetry of the KPM vectors used in the expansion, and it has no effect that the superlattice is skewed, since the periodic boundaries are take periodic in k, just like the vectors.

params["disorder"] = 0

kpm_params, energies = init_kpm_params(fsyst, params)

tile = kwant.TranslationalSymmetry(*honeycomb.prim_vecs)

%%time
with ProcessPoolExecutor() as executor:
    densities_unfolded = np.array(list(executor.map(
        density_in_k,
        k_path,
        [tile] * num_k_full,
    )))

CPU times: user 142 ms, sys: 53.7 ms, total: 196 ms
Wall time: 29.5 s

%%time
with ProcessPoolExecutor() as executor:
    densities_folded = np.array(list(executor.map(
        density_in_k,
        k_path,
        [fsyst._wrapped_symmetry] * num_k_full,
    )))

CPU times: user 138 ms, sys: 36.7 ms, total: 174 ms
Wall time: 21.7 s

plot_kpath_unfolded = hv.QuadMesh(
    (k_path_lengths, energies, densities_unfolded.T.real),
    kdims=[r'$k$', r'$\epsilon$'],
    vdims=hv.Dimension(r'$\rho(\epsilon, k)$', range=(0, 50)),
).opts(
    xticks=xticks,
    colorbar=True,
    cmap='bone_r',
    title="k-path unfolded",
)

plot_kpath_folded = hv.QuadMesh(
    (k_path_lengths, energies, densities_folded.T.real),
    kdims=[r'$k$', r'$\epsilon$'],
    vdims=hv.Dimension(r'$\rho(\epsilon, k)$', range=(0, 50)),
).opts(
    xticks=xticks,
    colorbar=True,
    cmap='bone_r',
    title="k-path folded",
)

plot_kpath_folded + plot_kpath_unfolded

Zone-unfolded k-path plot

Now we can add noise to the system. This will break the translational symmetry of the lattice, and only the translational symmetry of the superlattice will remain unbroken.

In this case, we expect to see some avoided crossing evidencing the zone-unfolding procedure.

params["disorder"] = 1

kpm_params, energies = init_kpm_params(fsyst, params)

tile = kwant.TranslationalSymmetry(*honeycomb.prim_vecs)

%%time
with ProcessPoolExecutor() as executor:
    densities_unfolded = np.array(list(executor.map(
        density_in_k,
        k_path,
        [tile] * num_k_full,
    )))

CPU times: user 134 ms, sys: 104 ms, total: 238 ms
Wall time: 25.5 s

%%time
with ProcessPoolExecutor() as executor:
    densities_folded = np.array(list(executor.map(
        density_in_k,
        k_path,
        [fsyst._wrapped_symmetry] * num_k_full,
    )))

CPU times: user 147 ms, sys: 35.5 ms, total: 183 ms
Wall time: 21 s

plot_kpath_unfolded = hv.QuadMesh(
    (k_path_lengths, energies, densities_unfolded.T.real),
    kdims=[r'$k$', r'$\epsilon$'],
    vdims=hv.Dimension(r'$\rho(\epsilon, k)$', range=(0, 50)),
).opts(
    xticks=xticks,
    colorbar=True,
    cmap='bone_r',
    title="k-path unfolded",
)

plot_kpath_folded = hv.QuadMesh(
    (k_path_lengths, energies, densities_folded.T.real),
    kdims=[r'$k$', r'$\epsilon$'],
    vdims=hv.Dimension(r'$\rho(\epsilon, k)$', range=(0, 50)),
).opts(
    xticks=xticks,
    colorbar=True,
    cmap='bone_r',
    title="k-path folded",
)

plot_kpath_folded + plot_kpath_unfolded