ARPES Simulation

In addition to the core ARPES_lib library, several other scripts in the module are written with the express purpose of facilitating calculation of the ARPES intensity. All relevant docs are included below.

In setting up an ARPES calculation, one requires an existing model system, i.e. instance of the TB_lib.TB_model class. This includes all relevant information regarding the orbital basis and the model Hamiltonian. In addition to this, a number of experimental parameters should be specified. Similar to the input for defining an orbital basis and a Hamiltonian, we use python dictionaries to specify these details. An example input is shown here.

ARPES_dict = {'cube':{'X':[-1,1,200],'Y':[-0.5,0.5,100],'E':[-1,0.05,1000],'kz':0},
                'hv':21.2,
                'pol':np.array([1,0,0]),
                'T':4.2,
                'SE':['fixed',0.02],
                'resolution':{'dE':0.005,'dk':0.01}}

By passing this dictionary to the input statement

experiment = ARPES_lib.experiment(TB,ARPES_dict)

We initialize an ARPES experiment which will diagonalize the Hamiltonian over a 200x100 point mesh in the region \(-1 \leq k_x \leq 1\) and \(-0.5 \leq k_y\leq0.5\) . For states in the range of energy \(-1\leq E\leq 0.05\) eV of the Fermi level, we will explicitly compute the ARPES cross section when requested. The calculation will be done with photon energy for He-1 \(\alpha\) at 21.2 eV, and a sample temperature of 4.2 K. This carries into the calculation in the form of a Fermi-Dirac distribution, which suppresses intensity from states with positive energies. The polarization of light is taken to be along the x direction, and is indicated by length-3 array associated with the keyword pol.

The SE key indicates the form of the self-energy to be imposed in evaluating the lineshape of the spectral features associated with each peak. Here we impose a fixed-linewidth of 20 meV on all states, to allow us to focus on the matrix elements alone, without further complications from energy-dependent lineshape broadening. As detailed in the ARPES_lib.SE_gen() below, more sophisticated options are available.

Finally, resolution is passed as well, with arguments for both energy and momentum resolution, expressed as full-width half maximum values, in units of electron volts and inverse Angstrom. Many more non-default options are available, including sample rotations, spin-projection (for spin-ARPES) and radial-integrals. See the documentation for ARPES_lib.experiment for further details.

Once a calculation of the matrix elements is completed, one is interested in plotting the associated intensity map. There are several options for this. First, the intensity map must be built using ARPES_lib.experiment.spectral(). Note that GUI tools do this automatically, without the user’s input. The ARPES_lib.experiment.spectral() method will apply the matrix element associated with each state, to its spectral function, and sum over all states to produce a complete dataset. It generates a raw and resolution broadened intensity map. The slice_select option allows for plotting specific cuts in energy or momentum. Once a full dataset has been generated, this can be passed to ARPES_lib.experiment.plot_intensity_map() to quickly image a different cut. Alternatively, interactive GUI tools are available under Matplotlib Plotter, or if the user has Tkinter installed, ARPES_lib.experiment.plot_gui().

Numerical Integration

For evaluation of radial integrals

\[B_{n,l}^{l'}(k) = (i)^{l'}\int dr R_{n,l}(r) r^3 j_{l'}(kr)\]

we use an adaptive integration algorithm which allows for precise and accurate evaluation of numeric integrals, regardless of local curvature. We do this by defining a partition of the integration domain which is recursively refined to sample regions of high curvature more densely. This is done until the integral converges to within a numerical tolerance.

The user is given some opportunity to specify details of the evaluation of radial integrals \(B_{n,l}^{l'}(k)\) used in the calculations. Specifications can be passed to the ARPES calculation through the ARPES_dict argument passed to the ARPES_lib.experiment object.

adaptive_int.general_Bnl_integrand(func, kn, lp)

Standard form of executable integrand in the e.r approximation of the matrix element

args:

  • func: executable function of position (float), in units of Angstrom
  • kn: float, norm of the k vector (in inverse Angstrom)
  • lp: int, final state angular momentum quantum number

return:

  • executable function of float (position)

***

adaptive_int.integrate(func, a, b, tol)

Evaluate the integral of func over the domain covered by a, b. This begins by seeding the evaluation with a maximally coarse approximation to the integral.

args:

  • func: executable
  • a: float, start of interval
  • b: float, end of interval
  • tol: float, tolerance for convergence

return:

  • Q: (complex) float, value of the integral

***

adaptive_int.rect(func, a, b)

Approximation to contribution of a finite domain to the integral, evaluated as a rough rectangle

args:

  • func: executable to evaluate
  • a: float, start of interval
  • b: float, end of interval

return:

  • recsum: (complex) float approximated area of the region under

function between a and b

***

adaptive_int.recursion(func, a, b, tol, currsum)

Recursive integration algorithm–rect is used to approximate the integral under each half of the domain, with the domain further divided until result has converged

args:

  • func: executable
  • a: float, start of interval
  • b: float, end of interval
  • tol: float, tolerance for convergence
  • currsum: (complex) float, current evaluation for the integral

return:

  • recursive call to the function if not converged, otherwise the result as complex (or real) float

***

radint_lib.define_radial_wavefunctions(rad_dict, basis)

Define the executable radial wavefunctions for computation of the radial integrals

args:

  • rad_dict: essential key is ‘rad_type’, if not passed,

assume Slater orbitals.

  • rad_dict[‘rad_type’]:

    • ‘slater’: default value, if ‘rad_type’ is not passed,
    Slater type orbitals assumed and evaluated for the integral
    • ‘rad_args’: dictionary of float, supplying optional final-state

    phase shifts, accounting for scattering-type final states. keys of form ‘a-n-l-lp’. Radial integrals will be accordingly multiplied

    • ‘hydrogenic’: similar in execution to ‘slater’,
    but uses Hydrogenic orbitals–more realistic for light-atoms
    • ‘rad_args’: dictionary of float, supplying optional final-state

    phase shifts, accounting for scattering-type final states. keys of form ‘a-n-l-lp’. Radial integrals will be accordingly multiplied

    • ‘grid’: radial wavefunctions evaluated on a grid of

    radii. Requires also another key_value pair:

    • ‘rad_args’: dictionary of numpy arrays evaluating

    the radial wavefunctions. Requires an ‘r’ array, as well as ‘a-n-l’ indicating ‘atom-principal quantum number-orbital angular momentum’. Must pass such a grid for each orbital in the basis!

    • ‘exec’: executable functions for each ‘a-n-l’ i.e.

    ‘atom-principal quantum number-orbital angular momentum’. If executable is chosen, require also:

    • ‘rad_args’, which will be a dictionary of

    executables, labelled by the keys ‘a-n-l’. These will be passed to the integral routine. Note that it is required that the executables are localized, i.e. vanishing for large radii.

    • ‘fixed: radial integrals taken to be constant float,

    require dictionary:

    • ‘rad_args’ with keys ‘a-n-l-lp’, i.e.

    ‘atom-principal quantum number-orbital angular momentum-final state angular momentum’ and complex float values for the radial integrals.

  • basis: list of orbital objects

return:

orbital_funcs: dictionary of executables

***

radint_lib.fill_radint_dic(Eb, orbital_funcs, hv, W=0.0, phase_shifts=None, fixed=False)

Function for computing dictionary of radial integrals. Can pass either an array of binding energies or a single binding energy as a float. In either case, returns a dictionary however the difference being that the key value pairs will have a value which is itself either a float, or an interpolation mesh over the range of the binding energy array. The output can then be used by either writing Bdic[‘key’] or **Bdic[‘key’]**(valid float between endpoints of input array)

args:

  • Eb: float or tuple indicating the extremal energies
  • orbital_funcs: dictionary of executable orbital radial wavefunctions
  • fixed: bool, if True, constant radial integral for each scattering

channel available: then the orbital_funcs dictionary already has the radial integral evaluated

  • hv: float, photon energy of incident light.

kwargs:

  • W: float, work function
  • phase_shifts: dictionary for final state phase shifts, as an optional

extension beyond pure- free electron final states. For now, float type.

return:

  • Brad: dictionary of executable interpolation grids

***

radint_lib.find_cutoff(func)

Find a suitable cutoff lengthscale for the radial integration: Evaluate the function over a range of 20 Angstrom, with reasonable detail (dr = 0.02 A). Find the maximum in this range. The cutoff tolerance is set to 1/1e4 of the maximum value. Since this ‘max’ is actually a lower bound on the true maximum, this will only give us a more strict cutoff tolerance than is absolutely possible. With this point found, we then find all points which are within the tolerance of zero. The frequency of these points is then found. When the frequency is constant and 1 for all subsequent points, we have found the point of convergence. If the ‘point of convergence’ is the last point in the array, the radial wavefunction really isn’t suitably localized and the user should not proceed without giving more consideration to the application of the LCAO approximation to such a function.

args:

  • func: the integrand executable

return:

  • float, cutoff distance for integration

***

radint_lib.gen_const(val)

Create executable function returning a constant value

args:

  • val: constant value to return when executable function

return:

  • lambda function with constant value

***

radint_lib.gen_orb_labels(basis)

Simple utility function for generating a dictionary of atom-n-l:[Z, orbital label] pairs, to establish which radial integrals need be computed.

args:

  • basis: list of orbitals in basis

return:

  • orbitals: dictionary of radial integral pairs

***

radint_lib.make_radint_pointer(rad_dict, basis, Eb)

Define executable radial integral functions, and store in a pointer-integer referenced array. This allows for fewer executions of the interpolation function in the event where several orbitals in the basis share the same a,n,l. Each of these gets 2 functions for l +/-1, which are stored in the rows of the array B_array. The orbitals in the basis then are matched to these executables, with the corresponding executable row index saved in B_pointers.

Begin by defining the executable radial wavefunctions, then perform integration at several binding energies, finally returning an interpolation of these integrations.

args:

  • rad_dict: dictionary of ARPES parameters: relevant keys are

‘hv’ (photon energy), ‘W’ (work function), and the rad_type (radial wavefunction type, as well as any relevant additional pars, c.f. radint_lib.define_radial_wavefunctions). Note: ‘rad_type’ is optional, (as is rad_args, depending on choice of radial wavefunction.)

  • basis: list of orbitals in the basis
  • Eb: tuple of 2 floats indicating the range of energy of

interest (increasing order)

return:

  • B_array: numpy array of Nx2 executable functions of float
  • B_pointers: numpy array of integer indices matching orbital

basis ordering to the functions in B_array

***

radint_lib.radint_calc(k_norm, orbital_funcs, phase_shifts=None)

Compute dictionary of radial integrals evaluated at a single |k| value for the whole basis. Will avoid redundant integrations by checking for the presence of an identical dictionary key. The integration is done as a simple adaptive integration algorithm, defined in the adaptive_int library.

args:

  • k_norm: float, length of the k-vector

(as an argument for the spherical Bessel Function)

  • orbital_funcs: dictionary, radial wavefunction executables

kwargs:

  • phase_shifts: dictionary of phase shifts, to convey final state scattering

returns:

  • Bdic: dictionary, key value pairs in form – ‘ATOM-N-L’:Bval

***

radint_lib.radint_dict_to_arr(Bdict, basis)

Take a dictionary of executables defined for different combinations of a,n,l and send them to an array, with a corresponding pointer array which can be used to dereference the relevant executable.

args:

  • Bdict: dictionary of executables with ‘a-n-l-l’’ keys
  • basis: list of orbital objects

return:

  • Blist: numpy array of the executables, organized by a-n-l,

and l’ (size Nx2, where N is the length of the set of distinct a-n-l triplets)

  • pointers: numpy array of length (basis), integer datatype

indicating the related positions in the Blist array

***

ARPES Library

ARPES_lib.G_dic()

Initialize the gaunt coefficients associated with all possible transitions relevant

return:
  • Gdict: dictionary with keys as a string representing (l,l’,m,dm) “ll’mdm” and values complex float.

All unacceptable transitions set to zero.

ARPES_lib.Gmat_make(lm, Gdictionary)

Use the dictionary of relevant Gaunt coefficients to generate a small 2x3 array of float which carries the relevant Gaunt coefficients for a given initial state.

args:
  • lm: tuple of 2 int, initial state orbital angular momentum and azimuthal angular momentum
  • Gdictionary: pre-calculated dictionary of Gaunt coefficients, with key-values associated with “ll’mdm”
return:
  • mats: numpy array of float 2x3
ARPES_lib.all_Y(basis)

Build L-M argument array input arguments for every combination of l,m in the basis. The idea is for a given k-point to have a single call to evaluate all spherical harmonics at once. The pointer array orb_point is a list of lists, where for each projection in the basis, the integer in the list indicates which row (first axis) of the Ylm array should be taken. This allows for very quick access to the l+/-1, m+/-1,0 Ylm evaluation required.

args:
  • basis: list of orbital objects
return:
  • l_args: numpy array of int, of shape len(lm_inds),3,2, with the latter two indicating the final state orbital angular momentum
  • m_args: numpy array of int, of shape len(lm_inds),3,2, with the latter two indicating the final state azimuthal angular momentum
  • g_arr: numpy array of float, shape len(lm_inds),3,2, providing the related Gaunt coefficients.
  • orb_point: numpy array of int, matching the related sub-array of l_args, m_args, g_arr related to each orbital in basis
ARPES_lib.con_ferm(ekbt)

Typical values in the relevant domain for execution of the Fermi distribution will result in an overflow associated with 64-bit float. To circumvent, set fermi-function to zero when the argument of the exponential in the denominator is too large.

args:
  • ekbt: float, (E-u)/kbT in terms of eV
return:
  • fermi: float, evaluation of Fermi function.
class ARPES_lib.experiment(TB, ARPES_dict)

The experiment object is at the centre of the ARPES matrix element calculation.This object keeps track of the experimental geometry as well as a local copy of the tight-binding model and its dependents. Such a copy is used to avoid corruption of these objects in the global space during a given run of the ARPES experiment.

args:

  • TB: instance of a tight-binding model object

  • ARPES_dict: dictionary of relevant experimental parameters including

    • ‘hv’: float, photon energy (eV),

    • ‘mfp’: float, mean-free path (Angstrom),

    • ‘resolution’: dictionary for energy and momentum resolution:

      • ‘dE’: float, energy resolution (FWHM eV),
      • ‘dk’: float, momentum resolution (FWHM 1/Angstrom)

    • ‘T’: float, Temperature of sample (Kelvin)

    • ‘cube’: dictionary momentum and energy domain

    (‘kz’ as float, all others ( ‘X’ , ‘Y’ , ‘E’ ) are list or tuple of floats Xo,Xf,dX)

optional args:

In addition to the keys above, ARPES_dict can also be fed the following:

  • ‘spin’: spin-ARPES measurement, list [+/-1,np.array([a,b,c])]

with the numpy array indicating the spin-projection direction (with respect to) experimental frame.

  • ‘rad_type’: string, radial wavefunctions, c.f. chinook.rad_int.py for details
  • ‘threads’: int, number of threads on which to calculate the matrix elements.

Requires very large calculation to see improvement over single core.

  • ‘slab’: boolean, will truncate the eigenfunctions beyond the penetration depth (specifically 4x penetration depth), default is False
  • ‘angle’: float, rotation of sample about normal emission i.e. z-axis (radian), default is 0.0
  • ‘W’: float, work function (eV), default is 4.0

***

M_compute(i)

The core method called during matrix element computation.

args:
  • i: integer, index and energy of state
return:
  • Mtmp: numpy array (2x3) of complex float corresponding to the matrix element projection for dm = -1,0,1 (columns) and spin down or up (rows) for a given state in k and energy.
Mk_wrapper(ilist)

Wrapper function for use in multiprocessing, to run each of the processes as a serial matrix element calculation over a sublist of state indices.

args:
  • ilist: list of int, all state indices for execution.
return:
  • Mk_out: numpy array of complex float with shape (len(ilist), 2,3)
SE_gen()

Self energy arguments are passed as a list, which supports mixed-datatype. The first entry in list is a string, indicating the type of self-energy, and the remaining entries are the self-energy.

args:
  • SE_args: list, first entry can be ‘func’, ‘poly’, ‘constant’, or ‘grid’

indicating an executable function, polynomial factors, constant, or a grid of values

return:
  • SE, numpy array of complex float, with either shape of the datacube,

or as a one dimensional array over energy only.

T_distribution()

Compute the Fermi-distribution for a fixed temperature, over the domain of energy of interest

return:
  • fermi: numpy array of float, same length as energy domain array defined by cube[2] attribute.
datacube(ARPES_dict=None, diagonalize=False)

This function computes the photoemission matrix elements. Given a kmesh to calculate the photoemission over, the mesh is reshaped to an nx3 array and the Hamiltonian diagonalized over this set of k points. The matrix elements are then calculated for each of these E-k points

kwargs:
  • ARPES_dict: can optionally pass a dictionary of experimental parameters, to update those defined

in the initialization of the experiment object.

return:
  • boolean, True if function finishes successfully.
diagonalize(diagonalize)

Diagonalize the Hamiltonian over the desired range of momentum, reshaping the band-energies into a 1-dimensional array. If the user has not selected a energy grain for calculation, automatically calculate this.

return:
None, however experiment attributes X, Y, ph, TB.Kobj, Eb, Ev, cube are modified.
gen_all_pol()

Rotate polarization vector, as it appears for each angle in the experiment. Assume that it only rotates with THETA_y (vertical cryostat), and that the polarization vector defined by the user relates to centre of THETA_x axis. Right now only handles zero vertical rotation (just tilt)

return:
  • numpy array of len(expmt.cube[1]) x 3 complex float, rotated polarization vectors

expressed in basis of spherical harmonics

plot_gui()

Generate the Tkinter gui for exploring the experimental parameter-space associated with the present experiment.

args:
  • ARPES_dict: dictionary of experimental parameters, c.f. the

__init__ function for details.

return:
  • Tk_win: Tkinter window.
plot_intensity_map(plot_map, slice_select, plot_bands=False, ax_img=None, colourmap=None)

Plot a slice of the intensity map computed in spectral. The user selects either an array index along one of the axes, or the fixed value of interest, allowing either integer, or float selection.

args:
  • plot_map: numpy array of shape (self.cube[0],self.cube[1],self.cube[2]) of float
  • slice_select: list of either [int,int] or [str,float], corresponding to

dimension, index or label, value. The former option takes dimensions 0,1,2 while the latter can handle ‘x’, ‘kx’, ‘y’, ‘ky’, ‘energy’, ‘w’, or ‘e’, and is not case-sensitive.

  • plot_bands: boolean, option to overlay a constant-momentum cut with

the dispersion calculated from tight-binding

  • ax_img: matplotlib Axes, for option to plot onto existing Axes
  • colourmap: matplotlib colourmap

return:

  • ax_img: matplotlib axis object
rot_basis()

Rotate the basis orbitals and their positions in the lab frame to be consistent with the experimental geometry

return:
  • list of orbital objects, representing a rotated version of the original basis if the

angle is finite. Otherwise, just return the original basis.

sarpes_projector()

For use in spin-resolved ARPES experiments, project the computed matrix element values onto the desired spin-projection direction. In the event that the spin projection direction is not along the standard out-of-plane quantization axis, we rotate the matrix elements computed into the desired basis direction.

return:
  • spin_projected_Mk: numpy array of complex float with same

shape as Mk

serial_Mk(indices)

Run matrix element on a single thread, directly modifies the Mk attribute.

args:
  • indices: list of all state indices for execution; restricting states
in cube_indx to those within the desired window
smat_gen(svector=None)

Define the spin-projection matrix related to a spin-resolved ARPES experiment.

return:
  • Smat: numpy array of 2x2 complex float corresponding to Pauli operator along the desired direction
spectral(ARPES_dict=None, slice_select=None, add_map=False, plot_bands=False, ax=None, colourmap=None)

Take the matrix elements and build a simulated ARPES spectrum. The user has several options here for the self-energy to be used, c.f. SE_gen() for details. Gaussian resolution broadening is the last operation performed, to be consistent with the practical experiment. slice_select instructs the method to also produce a plot of the designated slice through momentum or energy. If this is done, the function also returns the associated matplotlib.Axes object for further manipulation of the plot window.

kwargs:
  • ARPES_dict: dictionary, experimental configuration. See experiment.__init__ and experiment.update_pars()
  • slice_select: tuple, of either (int,int) or (str,float) format. If (int,int), first is axis index (0,1,2 for x,y,E) and the second is the index of the array. More useful typically is (str,float) format, with str as ‘x’, ‘kx’, ‘y’, ‘ky’, ‘E’, ‘w’ and the float the value requested. It will find the index along this direction closest to the request. Note the strings are not case-sensitive.
  • add_map: boolean, add intensity map to list of intensity maps. If true, a list of intensity objects is appended, otherwise, the intensity map is overwritten
  • plot_bands: boolean, plot bandstructure from tight-binding over the intensity map
  • ax: matplotlib Axes, only relevant if slice_select, option to pass existing Axes to plot onto
return:
  • I: numpy array of float, raw intensity map.
  • Ig: numpy array of float, resolution-broadened intensity map.
  • ax: matplotlib Axes, for further modifications to plot only if slice_select True
thread_Mk(N, indices)

Run matrix element on N threads using multiprocess functions, directly modifies the Mk attribute.

NOTE 21/2/2019 – this has not been optimized to show any measureable improvement over serial execution. May require a more clever way to do this to get a proper speedup.

args:
  • N: int, number of threads
  • indices: list of int, all state indices for execution; restricting

states in cube_indx to those within the desired window.

truncate_model()

For slab calculations, the number of basis states becomes a significant memory load, as well as a time bottleneck. In reality, an ARPES calculation only needs the small number of basis states near the surface. Then for slab-calculations, we can truncate the basis and eigenvectors used in the calculation to dramatically improve our capacity to perform such calculations. We keep all eigenvectors, but retain only the projection of the basis states within 2*the mean free path of the surface. The states associated with this projection are retained, while remainders are not.

return:
  • tmp_basis: list, truncated subset of the basis’ orbital objects
  • Evec: numpy array of complex float corresponding to the truncated eigenvector
array containing only the surface-projected wavefunctions
update_pars(ARPES_dict, datacube=False)

Several experimental parameters can be updated without re-calculating the ARPES intensity explicitly. Specifically here, we can update resolution in both energy and momentum, as well as temperature, spin-projection, self-energy function, and polarization.

args:

  • ARPES_dict: dictionary, specifically containing

    • ‘resolution’: dictionary with ‘E’:float and ‘k’:float
    • ‘T’: float, temperature, a negative value will suppress the Fermi function
    • ‘spin’: list of [int, numpy array of 3 float] indicating projection and spin vector
    • ‘SE’: various types accepted, see SE_gen for details
    • ‘pol’: numpy array of 3 complex float, polarization of light
kwargs:
  • datacube: bool, if updating in spectral, only the above can be changed. If instead, updating
at the start of datacube, can also pass:
  • hv: float, photon energy, eV
  • ang: float, sample orientation around normal, radiants
  • rad_type: string, radial integral type
  • rad_args: various datatype, see radint_lib for details
  • kz: float, out-of-plane momentum, inverse Angstrom
  • mfp: float, mean-free path, Angstrom
write_Ik(filename, mat)

Function for producing the textfiles associated with a 2 dimensional numpy array of float

args:

  • filename: string indicating destination of file
  • mat: numpy array of float, two dimensional
return:
  • boolean, True
write_map(_map, directory)

Write the intensity maps to a series of text files in the indicated directory.

args:
  • _map: numpy array of float to write to file
  • directory: string, name of directory + the file-lead name
return:
  • boolean, True
write_params(Adict, parfile)

Generate metadata text file associated with the saved map.

args:
  • Adict: dictionary, ARPES_dict same as in above functions, containing

relevant experimental parameters for use in saving the metadata associated with the related calculation.

  • parfile: string, destination for the metadata
ARPES_lib.find_mean_dE(Eb)

Find the average spacing between adjacent points along the dispersion calculated.

args:
  • Eb: numpy array of float, eigenvalues
return:
  • dE_mean: float, average difference between consecutive eigenvalues.
ARPES_lib.gen_SE_KK(w, SE_args)

The total self-energy is computed using Kramers’ Kronig relations:

The user can pass the self-energy in the form of either a callable function, a list of polynomial coefficients, or as a numpy array with shape Nx2 (with the first column an array of frequency values, and the second the values of a function). For the latter option, the user is responsible for ensuring that the function goes to zero at the tails of the domain. In the former two cases, the ‘cut’ parameter is used to impose an exponential cutoff near the edge of the domain to ensure this is the case. In all cases the input imaginary self-energy must be single-signed to ensure it is purely even function. It is forced to be negative in all cases to give a positive spectral function. With the input defined, along with the energy range of interest to the calculation, a MUCH larger domain (100x in the maximal extent of the energy region of interest) is defined wf. This is the domain over which we evaluate the Hilbert transform, which itself is carried out using: the scipy.signal.hilbert() function. This function acting on an array f: H(f(x)) -> f(x) + i Hf(x). It relies on the FFT performed on the product of the sgn(w) and F(w) functions, and then IFFT back so that we can use this to extract the real part of the self energy, given only the input. args:

w – numpy array energy values for the spectral peaks used in the ARPES simulation SE_args – dictionary containing the ‘imfunc’ key value pair (values being either callable, list of polynomial prefactors (increasing order) or numpy array of energy and Im(SE) values)

– for the first two options, a ‘cut’ key value pair is also required to force the function to vanish at the boundary of the Hilbert transform integration window.

return: self energy as a numpy array of complex float. The indexing matches that of w, the spectral features to be plotted in the matrix element simulation.

ARPES_lib.pol_2_sph(pol)

return polarization vector in spherical harmonics – order being Y_11, Y_10, Y_1-1. If an array of polarization vectors is passed, use the einsum function to broadcast over all vectors.

args:
  • pol: numpy array of 3 complex float, polarization vector in Cartesian coordinates (x,y,z)
return:
  • numpy array of 3 complex float, transformed polarization vector.
ARPES_lib.poly(input_x, poly_args)

Recursive polynomial function.

args:
  • input_x: float, int or numpy array of numeric type, input value(s) at which to evaluate the polynomial
  • poly_args: list of coefficients, in INCREASING polynomial order i.e. [a_0,a_1,a_2] for y = a_0 + a_1 * x + a_2 *x **2
return:
  • recursive call to poly, if poly_args is reduced to a single value, return explicit evaluation of the function.

Same datatype as input, with int changed to float if poly_args are float, polynomial evaluated over domain of input_x

ARPES_lib.progress_bar(N, Nmax)

Utility function, generate string to print matrix element calculation progress.

args:
  • N: int, number of iterations complete
  • Nmax: int, total number of iterations to complete
return:
  • st: string, progress status
ARPES_lib.projection_map(basis)

In order to improve efficiency, an array of orbital projections is generated, carrying all and each orbital projection for the elements of the model basis. As these do not in general have the same length, the second dimension of this array corresponds to the largest of the sets of projections associated with a given orbital. This will in practice remain a modest number of order 1, since at worst we assume f-orbitals, in which case the projection can be no larger than 7 long. So output will be at worst len(basis)x7 complex float

args:
  • basis: list of orbital objects
return:
  • projarr: numpy array of complex float

Intensity Maps

class intensity_map.intensity_map(index, Imat, cube, kz, T, hv, pol, dE, dk, self_energy=None, spin=None, rot=0.0, notes=None)

Class for organization and saving of data, as well as metadata related to a specific ARPES calculation.

copy()

Copy-by-value of the intensity map object.

return:

  • intensity_map object with identical attributes to self.

***

save_map(directory)

Save the intensity map: if 2D, just a single file, if 3D, each constant-energy slice is saved separately. Saved as .txt file

args:

  • directory: string, directory for saving intensity map to

return:

  • boolean

***

write_2D_Imat(filename, index)

Sub-function for producing the textfiles associated with a 2dimensional numpy array of float

args:

  • filename: string, indicating destination of file
  • index: int, energy index of map to save, if -1, then just a 2D map, and save the whole

thing

return:

  • boolean

***

write_meta(destination)

Write meta-data file for ARPES intensity calculation.

args:

  • destination: string, file-lead

return:

  • boolean

***

Matplotlib Plotter

A built-in data-explorer is included in chinook, built using matplotlib to ensure cross platform stability. The figure below shows an example screen capture for a calculation on \(Sr_2IrO_4\) . The user has the ability to scan through the momentum and energy axes of the dataset, and the cursor can be used actively to select momentum- and energy- distribution curves in the side and lower panels. A scatterplot of the bare dispersion, as computed from the Hamiltonian diagonalization is plotted overtop the intensity map.

_images/matplotlib_plotter.png

Created on Thu Sep 12 14:53:36 2019

@author: rday

class matplotlib_plotter.interface(experiment)

This interactive tool is intended for exploring the dataset associated with an ARPES simulation using chinook. The user can scan through the datacube in each dimension. This uses matplotlib natively rather than alternative gui systems in python like Tkinter, which makes it a bit more robust across platforms.

bin_energy()

Translate the exact energy value for the band peaks into the discrete binning of the intensity map, to allow for cursor queries to be processed.

return:

  • coarse_pts: numpy array of float, same lengths as self.state_coords,

but sampled over a discrete grid.

***

find_cursor()

Find nearest point to the desired cursor position, as clicked by the user. The cursor event coordinates are compared against the peak positions from the tight-binding calculation, and a best choice within the plotted slice is selected.

args:

  • cursor: tuple of 2 float, indicating the column and row of the

event, in units of the data-set scaling (e.g. 1/A or eV)

***

plot_img()

Update the plotted intensity map slice. The plotted bandstructure states are also displayed.

***

run_gui()

Execution of the matplotlib gui. The figure is initialized, along with all widgets and chosen datasets. The user has access to both the slice of ARPES data plotted, in addition to the orbital projection plotted in upper right panel.

Tilt Geometry

One can account for rotation of the experimental geometry during acquisition of a particular dataset by performing an inverse rotation on the incoming light vector. Similarly, spin projections can also be rotated into the laboratory frame to reflect the effect of a misaligned or rotated sample in a spin-ARPES experiment.

Created on Fri Dec 28 13:33:59 2018

@author: ryanday

tilt.ang_mesh(N, th, ph)

Generate a mesh over the indicated range of theta and phi, with N elements along each of the two directions

args:

  • N: int or iterable of length 2 indicating # of points along th, and ph respectively
  • th: iterable length 2 of float (endpoints of theta range)
  • ph: iterable length 2 of float (endpoints of phi range)

return:

  • numpy array of N_th x N_ph float, representing mesh of angular coordinates

***

tilt.gen_kpoints(ek, N, thx, thy, kz)

Generate a mesh of kpoints over a mesh of emission angles.

args:

  • ek: float, kinetic energy, eV
  • N: tuple of 2 int, number of points along each axis
  • thx: tuple of 2 float, range of horizontal angles, radian
  • thy: tuple of 2 float, range of vertical angles, radian
  • kz: float, k-perpendicular of interest, inverse Angstrom

return:

  • k_array: numpy array of N[1]xN[0] float, corresponding to mesh of in-plane momenta

***

tilt.k_mesh(Tmesh, Pmesh, ek)

Application of rotation to a normal-emission vector (i.e. (0,0,1) vector) Third column of a rotation matrix formed by product of rotation about vertical (ky), and rotation around kx axis c.f. Labbook 28 December, 2018

args:

  • Tmesh: numpy array of float, output of ang_mesh
  • Pmesh: numpy array of float, output of ang_mesh
  • ek: float, kinetic energy in eV

return:

  • kvec: numpy array of float, in-plane momentum array associated with angular emission coordinates

***

tilt.k_parallel(ek)

Convert kinetic energy in eV to inverse Angstrom

tilt.plot_mesh(ek, N, th, ph)

Plotting tool, plot all points in mesh, for an array of N angles, at a fixed kinetic energy.

args:

  • ek: float, kinetic energy, eV
  • N: tuple of 2 int, number of points along each axis
  • thx: tuple of 2 float, range of horizontal angles, radian
  • thy: tuple of 2 float, range of vertical angles, radian

***

tilt.rot_vector(vector, th, ph)

Rotation of vector by theta and phi angles, about the global y-axis by theta, followed by a rotation about the LOCAL x axis by phi. This is analogous to the rotation of a cryostat with a vertical-rotation axis (theta), and a sample-mount tilt angle (phi). NOTE: need to extend to include cryostats where the horizontal rotation axis is fixed, as opposed to the vertical axis–I have never seen such a system but it is of course totally possible.

args:

  • vector: numpy array length 3 of float (vector to rotate)
  • th: float, or numpy array of float – vertical rotation angle(s)
  • ph: float, or numpy array of float – horizontal tilt angle(s)

return:

  • numpy array of float, rotated vectors for all angles: shape 3 x len(ph) x len(th)

NOTE: will flatten any length-one dimensions

***