Posts Matplotlib: 1D Chain ARPES Signal
Post
Cancel

Matplotlib: 1D Chain ARPES Signal

ARPES Signal of 1D Chain within Tight-binding 1

  • The eigenvalues of 1D mono-atomic chain within the tight-binding model

    \[\begin{equation*} \varepsilon_\kappa = \varepsilon_0 - 2t\cos(\frac{\kappa}{N+1}\pi)\qquad \kappa=1,\ldots,N, \end{equation*}\]
  • The ARPES signal of the linear chain from Eq.(54) of Ref1

    \[\begin{equation*} w_{fi} = \sum_\kappa \bigl| \langle \mathbf{k}_f | \mathbf{R},\kappa \rangle \bigr|^2 \, {\cal G}_\kappa(\omega - \varepsilon_\kappa) \end{equation*}\]

    where ${\cal G}_\kappa$ is a Gaussian broadening function and $\langle\mathbf{k}_f \vert\mathbf{R},\kappa\rangle$ is given by

    \[\begin{equation*} \langle \mathbf{k}_f | \mathbf{R},\kappa \rangle \propto \sum_j^N \sin\left( \frac{j\kappa}{N+1}\pi \right) \, e^{-\imath k_{fx}\cdot x_j} \end{equation*}\]

Code

Below is the code and figure that simulate the ARPES signals, similar to Fig. 9 in Ref 1.

Figure 1. ARPES signal of the 1D linear chain with different lengths. Figure generated by Matplotlib.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#!/usr/bin/env python

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator

plt.style.use('dark_background')
mpl.rcParams['axes.unicode_minus'] = False

def photoem_mat_elem(k0, N=3):
    '''
    '''

    k0    = np.asarray(k0)
    nkpts = k0.size
    # photoemission matrix elements
    pk = np.array([
        np.sin(m*j*np.pi / (N+1))*np.exp(-1j*k0*j)
        for m in range(1, N+1)
        for j in range(1, N+1)
    ]).reshape((N, N, nkpts)).sum(axis=1)

    ek = np.array([
        -2*np.cos(m*np.pi/(N+1))
        for m in range(1, N+1)
    ])

    return ek, np.abs(pk)**2

if __name__ == "__main__":
    nkpts = 400      # No. of k-points
    nedos = 500      # No. of points for DOS
    sigma = 0.05     # energy broadening

    # k-space grid
    k0     = np.linspace(-np.pi, np.pi, nkpts)
    # energy grid
    e0     = np.linspace(-2-sigma*5, 2+sigma*5, nedos)
    x0, y0 = np.meshgrid(k0, e0, indexing='ij')

    ############################################################
    fig = plt.figure(
        figsize=(7.2, 4.8),
        dpi=480,
        constrained_layout=True
    )

    # 3x3 subplots
    axes_array = np.arange(9, dtype=int).reshape((3, 3))

    axes = fig.subplot_mosaic(
        axes_array,
        empty_sentinel=-1,
        gridspec_kw=dict(
            # height_ratios= [1, 0.5, 0.5],
            # width_ratios=[2, 2],
            # hspace=0.05,
            # wspace=0.06,
        )
    )
    axes = [axes[ii] for ii in range(axes_array.max()+1)]
    ############################################################

    # chain_lengths = np.arange(9) + 1
    # -1 for infinite chain length
    chain_lengths = np.array([1, 2, 3, 4, 5, 6, 10, 20, -1])

    for ii in range(chain_lengths.size):
        N  = chain_lengths[ii]
        ax = axes[ii]

        # finite chain length
        if N > 0:
            ek, pk = photoem_mat_elem(k0, N)
            # smearing
            ss = np.dot(
                pk.T,
                (1 / np.sqrt(2*np.pi) / sigma) * np.exp(-(e0[None, :] - ek[:, None])**2 / 2 / sigma**2)
            )
        # infinite chain length
        else:
            ek = -2 * np.cos(k0)
            # smearing
            ss = (1 / np.sqrt(2*np.pi) / sigma) * np.exp(-(e0[None, :] - ek[:, None])**2 / 2 / sigma**2)

        ax.pcolormesh(x0, y0, ss, cmap='magma')
        ax.plot(k0, -2*np.cos(k0), ls='--', lw=0.4)

        ax.yaxis.set_minor_locator(AutoMinorLocator(n=2))

        ax.set_xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi])
        ax.set_xticklabels([
            r'-$\,\frac{\pi}{a}$',
            r'-$\,\frac{\pi}{2a}$',
            '0',
            r'$\frac{\pi}{2a}$',
            r'$\frac{\pi}{a}$'
        ])

        if ii > 5:
            ax.set_xlabel(r'$k$')

        if ii % 3 == 0:
            ax.set_ylabel(r'Energy ($t$)', labelpad=5)

        if N > 0:
            ax.text(0.05, 0.05,
            # ax.text(0.50, 0.95,
                r'$N={}$'.format(N),
                # ha='center', va='top',
                ha='left', va='bottom',
                fontsize='small',
                transform=ax.transAxes
            )
        else:
            ax.text(0.05, 0.05,
            # ax.text(0.50, 0.95,
                r'$N=\infty$',
                # ha='center', va='top',
                ha='left', va='bottom',
                fontsize='small',
                transform=ax.transAxes
            )

    plt.savefig('arpes_1d_tb.png')
    from subprocess import call
    call('feh -xdF arpes_1d_tb.png'.split())

References

This post is licensed under CC BY 4.0 by the author.