Computational Framework for Angle-Resolved Photoemission Spectroscopy

We have developed the numerical software package $chinook$, designed for the simulation of photoemission matrix elements. This quantity encodes a depth of information regarding the orbital structure of the underlying wavefunctions from which photoemission occurs. Extraction of this information is often nontrivial, owing to the influence of the experimental geometry and photoelectron interference, precluding straightforward solutions. The $chinook$ code has been designed to simulate and predict the ARPES intensity measured for arbitrary experimental configuration, including photon-energy, polarization and spin-projection, as well as consideration of both surface-projected slab and bulk models. This framework then facilitates an efficient interpretation of the photoemission experiment, allowing for a deeper understanding of the electronic structure in addition to the design of new experiments which leverage the matrix element effects towards the objective of selective photoemission from states of particular interest.


INTRODUCTION
Angle-resolved photoemission spectroscopy (ARPES) and its variants have developed in recent years to be established among the pre-eminent experimental methods in solid-state physics. With an intimate connection to the one-electron removal spectral function, ARPES is unique among the suite of techniques available to condensed matter physicists in its direct correspondence to the electronic structure of crystalline materials, providing access to the one electron removal spectral function within its native momentum space [1][2][3].
In the framework of Fermi's Golden Rule, the photoemission intensity is described as: where A f,i is the one-electron removal spectral function, and | ψ f |∆ |ψ i | 2 the photoemission matrix element. The spectral function: (2) reflects the overlap between the initial N -particle manybody wavefunction upon removal of an electron and the ensemble of (N − 1)-particle final state wavefunctions, while preserving energy conservation. Written as the imaginary part of the retarded Green's function: The spectral function is seen to carry details of both the underlying bare dispersion associated with the electronic structure of the material 0 k , as well as correlations a The full chinook software package and documentation are available at https://www.github.com/rpday/chinook via the self energy Σ(k, ω) = Σ (k, ω) + iΣ"(k, ω). In the opposing limits of vanishing and strong interactions, ARPES is described as an ideal probe of the bandstructure and correlation effects respectively.
In practice, the photoemission can be strongly modulated by the | ψ f |∆ |ψ i | 2 term, altering the spectral intensity through the dependence of the initial and final states on energy, momentum, and band index. At worst, this suppresses all intensity from certain bands, precluding their study by ARPES entirely. From a different perspective however, this modulation can be viewed as an additional experimental signature in the ARPES intensity which encodes a description of the photoemitted electron's wavefunction.
This term can be simulated to allow for quantitative descriptions and insights regarding the experimental signal. While such an approach has been made at some level for a number of ARPES experiments, this requires substantial effort in developing a specific model for each study [4][5][6][7][8][9]. The development of a standard numerical framework would allow for a much larger set of experiments to be analysed at this level, providing the opportunity to understand and leverage the matrix element effects in a broad class of materials. We have pursued this objective through the development of an open source software package, chinook, implemented in Python to enable a broad audience to perform quick and easy simulation of photoemission-related phenomena, thereby improving both the interpretation and analysis of experimental data.
In the following, we will outline the primary workflow of our numerical approach, and the various ways in which this package can be applied to the study of the electronic structure of solids via ARPES.
FIG. 1. Experimental ARPES on FeSe. Both panels display ARPES intensity from valence states taken at hν = 37 eV and 120 K, directed along the ΓM direction. Polarization is set to linear vertical (a) and horizontal (b), allowing for photoemission from states of different orbital character, as indicated by the insets. Adapted from Ref. [11].

Matrix Element Effects
In the design of an ARPES experiment, a cursory understanding or prediction of the matrix elements relevant to a given system can dramatically improve one's ability to study aspects of the electronic structure of interest. Before proceeding to explicit description of our software, it is instructive to consider a motivating example, taken here to be the iron-based superconductor FeSe. In Fig.  1, we plot experimental data along the ΓM direction taken with two, orthogonal linear polarizations of light at hν = 37 eV. Near the Brillouin zone centre, three hole-bands disperse away from the Fermi level (c.f. Fig.  2). For light polarized along the momentum axis in Fig.  1a, only a single state is observed clearly, whereas the perpendicular polarization in Fig. 1b. illuminates this and several other states. The third hole band is almost imperceptible, for any choice of polarization. These observations can be explained through an understanding of the orbital structure of the underlying electronic states, indicated by the insets of Fig. 1, in combination with the orbital-mixing effects of spin-orbit coupling (SOC), as will be discussed in more detail below [10,11]. Recent experiments designed with these effects in mind have leveraged the dipole matrix elements to perform targeted spin-and angle-resolved photoemission from states of particular interest, extracting fundamental information pertaining to a broad variety of orbital-related phenomena [11][12][13][14][15][16].
In many cases, orbital symmetry can be extracted from the polarization dependence alone. The information available from these arguments is limited, in particular for high-orbital (l > 1) states, where the dimension of the vector field provides insufficient means to iden-tify all orbitals uniquely. This is further complicated in multi-atom bases, where now relative phases between different sites can differ from symmetric combinations in a momentum-dependent fashion. In these situations, the photoemission intensity pattern is found to depend sensitively on the relative phases within the initial state wavefunction, producing so-called photoelectron interference patterns. In this way however, the matrix elements encode further information regarding the initial state beyond orbital symmetry alone. These effects have been seen in for example graphene [7,17] and topological insulators [8,9], demonstrating the full depth of information regarding the initial state wavefunction which is contained in the ARPES matrix element. To further leverage the information available from ARPES experiments, it is advantageous to be able to simulate the full ARPES experimental intensity, while maintaining physical transparency. By preserving access to the relevant model parameters, one can then establish a more fundamental, and conceptual understanding of the electronic structure.

Model Hamiltonian
There are various levels at which the ARPES matrix element can be modeled [1,4,6,[18][19][20]. While the most sophisticated approaches account for the possibility of scattering processes subsequent to the photoemission event such as those which make use of Korringa-Kohn-Rostoker final states [21][22][23], we make two important simplifying assumptions here. First, the final states are taken to be free-electron plane waves: At high photon energies, the assumption of the plane wave final state is particularly well justified, as the crystal potential can be treated as a perturbation and sensitivity to the momentum structure of the exact final states becomes negligible [3]. The validity of this assumption is ultimately material dependent, however similar logic as that applied to the domain of suitability for the Born approximation can be made: such an assumption is reasonable when either the crystal potential V o 2 /m e a 2 or in the high-energy limit, ka, where a is the range of the potential. At present it is possible within chinook to relax this assumption only in the restricted sense of Ref. 18, as one can include phase shifts to the final state expansion. While beyond the scope of chinook in its current form, it would be possible to write the final states in the form of more sophisticated scattering final states, where the radial and orbital components of the ket in Eq. 4 are modified appropriately to reflect the presence of a finite crystal potential. This Fe-3d density of states is shown in (b). In (c) we plot the crystal structure projected into 2D, with Fe in red and the Se above (large) and below (small) in grey. In plane hopping terms are indicated, in addition to the primitive unit cell. The Fermi surface is plotted in (d).
can be done through modification of the radial integrals B l b defined below. Secondly, we work within a tight-binding framework wherein the initial states can be described by localized atomic-like orbitals, centred on the sites of the lattice basis. In materials where the spin degree of freedom is relevant, the orbital basis can be doubled to define a complete spinor basis, represented here by χ ± . Formally, the tight-binding basis set is expressed typically as: where a represents a basis index and n, l the principal and orbital quantum numbers respectively. R a n,l (r) is a radial wavefunction, K a l (Ω) a cubic harmonic, and χ ± the spinor projection. Alternatives such as distorted and rotated basis states can also be accommodated, so long as a unitary transformation into the basis of spherical harmonics can be made for the purpose of photoemission calculations. While these simplifications are in some cases unable to capture the full structure of the experimental photoemission intensity, we trade this level of universality for the substantial gains in transparency and physical insight which can be extracted from this approach.
Regarding the definition of the tight-binding model, there are various formalisms which are found in the literature, including Slater-Koster [25], t ab , and Wannier Hamiltonians: we have made an effort to accommodate all possible variations without loss of functionality. We require only that the Hamiltonian matrix elements can be written as a Fourier series of bilinear terms in the orbital Hilbert space: In this expression r ab denotes the full connecting vector between basis states φ a and φ b , as opposed to the equivalent form where one refers to the connecting lattice vector alone. In addition to H o (k), any other bilinear functions of momentum can also be added to the full Hamiltonian at this stage, including spin-orbit coupling, and orbital or spin order. Adherence to Bloch's theorem can be further relaxed in the context of the low-energy effective models which describe a narrow region of momentum space: in this scenario, a more general function of momentum which satisfies the point-group symmetry can be employed [26].
With the basis and Hamiltonian so defined, the eigenvalue problem can be solved, and the initial state wavefunctions then defined as a superposition of the basis states described by Eq. 5: With this information, a full characterization of the model for a specific material can be performed, followed by subsequent simulation of ARPES matrix elements. This includes density of states, bandstructure, 3D Fermi surface, orbital projections of the eigenstates, as well as the expectation values of various operators of interest, such as for example, L · S as in Fig. 2 for tetragonal FeSe. By defining an N ×N Hermitian matrix, the expectation value of any observable operator can be computed in this way.

Computation of Matrix Elements
The workflow of chinook is sketched in Fig. 3. Once a satisfactory material model is established, one can proceed to the simulation of ARPES intensity maps. A suitable region of interest in momentum and energy space must be defined, and the eigenvalue problem is then solved over this domain.
We model the matrix elements of the dipole operator as: where we have made use of the commutation relations to express the dipole operator in the position representation.
In explicit evaluation of the ARPES matrix element, we expand both the initial and final states as prescribed by Eqs. 4 and 5, which allows us to express Eq. 8 as: This sum over integrals can be expressed in compact form as: and µ the components of the polarization vector. In the third line, we have absorbed an extinction factor e −z b /2ξ into η b , where ξ represents the mean-free path of a photoemitted electron, and z b the spatial extent of the basis orbital below the surface [8,27].The radial and angular integrals are contained in the: and : terms respectively. Here, G b,µ l ,m is equivalent to a small subset of Gaunt coefficients, allowing for efficient and exact evaluation of this term. Meanwhile, the radial integrals can not necessarily be expressed in an analytical form and must be computed numerically, as the radial wavefunction is loosely constrained in most tight-binding models [28]. Whether a hydrogenic, Slater, or more complex object should be employed to describe the radial wavefunction is left to the discretion of the user, as the best choice is somewhat dependent on the nature of the material and states of interest. The user is given the opportunity to select from a variety of initial state wavefunctions in addition to importing their own functions or radial integrals at the start of the calculation. This could be the Wannier function grid as generated by for example Wannier90 [29].
We note that elsewhere it is common to take advantage of the plane-wave final state to recast the matrix element as a polarization modulated Fourier transform of the initial state [19,30]. Specifically, M F T can be expanded as in Eq. 10, establishing some formal equivalence to Eq. 8. However, one commonly observes qualitative deviation from experiment within this description, due to the form of the radial integrals B l b introduced above. In the Fourier representation, these are written as: One can contrast this with Eq. 11 which we employ in chinook. It is made explicit in Eq. 14 that the Fourier representation of the dipole operator imposes radial integrals which are independent of final state angular momentum. The implications for multi-orbital systems, and for those where l > 0, are important as final state interference becomes relevant. This is ultimately why for a plane-wave final state, the position, rather than momentum representation of the dipole operator yields a better description of experiment. We emphasize that although the limited constraints of tight-binding imply that the integrals B l b are to some extent parameters of the calculation, support for distinct final state angular momentum cross-sections is essential to the success of the position representation used here. Furthermore, it offers a natural extension of our framework to scattering-final states, wherein the commutation relations required to establish Eq. 8 are more rigourously justified. Further discussion of these approximations can be found in the Supplementary Materials.
Returning to the calculations executed in chinook, the central object of importance is the coherent matrix element factor: Evaluation of this object proceeds following Eq. 9. Each column corresponds to the projection of the polarization vector in the basis of spherical harmonics, (i.e. µ ≡ ∆m = ±1, 0), and the rows indicate the spinor projections. By retaining the matrix element in this coherent form, the ARPES intensity for arbitrary polarization and spin projection can be recalculated at run-time with minimal computational overhead. For each band and k-point in the region of interest, a spectral function as defined in Eq. 3 is added to the total intensity map, with its amplitude multiplied by: The photoemission intensity is then computed as described in Eq. 1, with | ψ f |∆ |ψ i | 2 → |M α | 2 . Spin projection, polarization, resolution, temperature and selfenergy can all be updated with little overhead at runtime.
With ARPES intensity maps then calculated for different experimental configurations, these results can be exported for further analysis, or combined to define quantities such as spin-polarization and circular/linear dichroism. In this sense, the output of the standard chinook calculation is a three-dimensional array of intensity in coordinates of momentum and energy which can be explored and analyzed in the same way as an experimental ARPES measurement.

APPLICATIONS Bulk Electronic Structure and Orbital Texture
Returning to the motivating case of the Fe-based superconductor FeSe, we can implement the model characterized above in Fig. 2 and compare the simulated ARPES intensity against the low energy region of Fig.  1. As with the experiments, the calculations were done at hν = 37 eV and T = 120 K. A Fermi-liquid type self-energy has been applied to the spectral features, resulting in an energy-dependent broadening of the photoemission linewidth. In the present case, the tight-binding model has already been renormalized to match the experimental spectra, such that the dispersion is more appropriately defined as k = 0 k − Σ (k, ω). Consequently, the self-energy used in the ARPES simulation is purely imaginary, Σ( k, ω) = iΣ"(ω) = −i(0.005 + 1.0ω 2 ), which is plotted in Fig. 4e. As ARPES matrix-elements can confound the evaluation of the spectral function and correlation effects in experimental data [31][32][33], the ability to model both components in the same environment can facilitate the disentanglement of these two objects of interest.
The simulation in Fig. 4 captures the relative intensity ratio between the three hole bands, with the heaviest (largest effective mass) band visible only through the SOC-induced hybridization gaps near E B = 50 meV. This latter state, composed primarily of d xy orbitals, has vanishing photoemission intensity along the normal emission (lim k →0 ) direction due to the selection rules associated with its definition in terms of spherical harmonics Y ±2 2 : all possible final states have a node along the normal emission direction. While conventional interpretation of the remaining states assumes d xz/yz -like wavefunctions, SOC allows for finite intensity from both states, as observed both experimentally and in the simulation near k = 0Å −1 . These inferences regarding the orbital structure of the initial states are supported by projection of the tight-binding eigenstates onto the basis of spherical harmonics, as done at select k-points in Fig.  4a using built-in diagnostic tools from chinook.
A more direct measure of the influence of SOC can be achieved through combining circularly polarized light with spin resolution to gain explicit access to both spin and orbital degrees of freedom. One can define the polarization asymmetry as: where subscripts indicate the helicity of light polarization, and superscripts the spin-projection of the photoelectrons. This quantity is closely related to the projection of L · S along the quantization axis of the experiment, allowing for a connection between Fig. 2c and Fig.  4d to be made. This technique has been applied to both ruthenates [13,34] and Fe-based superconductors [11], and utilizes the dipole selection rules encoded within the matrix element factor to provide the most direct measure of spin-orbital entanglement in solid state.

Supercell Impurity Model
Consideration of a supercell model illustrates the information encoded in the ARPES matrix element beyond orbital symmetry alone. Regarding the electronic structure of periodic systems, one can choose an arbitrarily large unit cell in exchange for a reduced Brillouin zone and additional backfolded bands. By contrast, impurities or other symmetry-breaking potentials (SBP) explicitly require such an expanded unit cell. While one can numerically perform an unfolding of these bands in an attempt to recover the spectrum within the extended Brillouin zone [35], such an unfolding is carried out naturally in the photoemission experiment.
In the absence of the SBP, Bloch's theorem would impose that the original band becomes a symmetric superposition over the neighbouring lattice sites. The additional bands, which must be orthogonal to the original state will destructively interfere in evaluation of the ARPES matrix element, preventing observation of many of the folded states. Ultimately, an SBP can mix these states; when the SBP is an essential feature of the potential landscape, as in graphene and the Fe-based superconductors [7,36], the folded bands can be observed with strong intensity over a range of momenta. When the SBP is weak or disordered, intensity from these folded bands vanishes away from the avoided crossings. To demonstrate these effects, we consider the artificial example of a square lattice of Li 2s orbitals, into which we substitute some number of Na 3s orbitals. Allowing for nearest neighbour hopping alone, and imposing a N a = -0.35 eV impurity potential for the Na sites, we simulate the effect of local defects in this lattice and the resulting ARPES spectra. Kinetic terms and onsite potentials have been adapted from the phenomenological rules set out in Ref. 37. In an attempt to consider the impurity problem realistically, we populate a 30×30 supercell of Li with various concentrations of randomly distributed Na impurities. For each particular distribution, the density of states is integrated to fix the Fermi-level at half-filling, consistent with electron counting. We then compute the photoemission intensity at hν = 21.2 eV over the ΓM direction of the extended Brillouin zone. For clarity, we assume a constant intrinsic linewidth of 10 meV along the entire dispersion. Energy and momentum resolution are set to 10 meV and 0.005Å −1 . The results, plotted in Fig. 5c have been averaged over 80 such configurations, corresponding to a nominal doping of Li 0.9 Na 0.1 . This can be compared against the pure Li-supercell in Fig. 5b. The full spectrum of the latter is displayed in Fig. 5a.
In each case, we compute the photoemission intensity from all states. However intensity from all folded bands is vanishing in the absense of the SBP; the full bandstructure is plotted in white over the spectra to demonstrate the large suppression of photoemission intensity. At the bottom of the band, the dipole selection rules suppress photoemission intensity from even the main band. In the disorder-averaged supercell, a substantial broadening of the spectral lineshape is observed [38,39]. As indicated by the overlain bandstructures of Fig. 5a-c, the impurity potential introduces a high density of avoided crossings, where the eigenvector supports finite photoemission intensity. In this sense, the broadening can be associated with the relative phases within the tight-binding eigenvector to which the ARPES matrix element is sensitive.
One can demonstrate that the linewidth broadening is dependent on both concentration and strength of impurities. In Fig. 5e,f we plot energy distribution curves (EDCs) at k F = 0.44Å −1 for fixed concentration (Li 0.9 Na 0.1 ) with variable attractive (negative) N a , and fixed N a = −0.35 eV with variable concentration. Each spectrum has been averaged over 80 similar configurations, and normalized to its peak intensity. The linewidth is observed to increase monotonically with both concentration and impurity potential, indicating the similar role these degrees of freedom play in modifying the spectra of disordered systems. Despite this lineshape broadening, the low-energy dispersion is resilient against a high level of disorder, as illustrated by Fig. 5c. By applying an out-of-plane polarization sensitive to states FIG. 6. Surface-projected photon energy dependence in FeSe. In (a), we build a 20 layer slab tight-binding model from the bulk model in Fig. 2. We calculate the ARPES intensity along the ΓM direction, at several photon energies; the results are summarized by a cut at constant binding energy of EB = 25 meV in panel (b). Spectra at each photon energy have been renormalized to their maximum intensity. In panels (c-e), spectra at select photon energies are plotted, chosen to correspond to the same kz = 0Å −1 point in successive Brillouin zones to enable direct comparison. Photon energy dependence of the measured linewidth is observed due to effective kz integration.
near the bottom of the band, we also confirm an increase of the bandwidth for this attractive impurity potential, which grows quadratically with impurity potential for the modest | N a | W < 0.15 considered here: at 10% Na and N a = −0.35 eV, the band bottom is extended 30 meV. Such detailed study of the impurity-substituted ARPES spectra is not possible without consideration of the matrix elements, which allow for a straightforward disentanglement of the supercell bandstructure and an opportunity to achieve meaningful insights from disordered materials.

Surface vs Bulk, and Emergent kz Dispersion
To this point we have considered the bulk-electronic structure, but it is important to appreciate the surfacesensitivity of the ARPES experiment: the high scattering cross-section in the ultraviolet regime results in penetration depths of the order of 5-10Å [27]. This corresponds to the top few unit cells of the lattice, depending on ex-perimental details. In many cases, the surface introduces modest corrections to the local electronic structure, facilitating a direct connection between the measured photoemission intensity and the bulk electronic structure [40]. In others, details of the surface preparation result in reconstructions of the ARPES spectra which deviate profoundly from the bulk electronic structure [40][41][42][43].
This surface sensitivity becomes rather important in the context of three-or even quasi-two-dimensional materials, where the photoelectron escape depth and k z information are intimately connected. For intermediate energies in the ultraviolet regime, where the penetration depth is of the order of 5Å, the ∆k z required by the uncertainty principle becomes comparable to the size of the Brillouin zone. In the presence of finite k z dispersion, this can result in anomalously broad linewidths, as the spectrum effectively integrates over the third dimension of momentum space. This is visualized well in Fig. 6, where we have projected our FeSe model onto a 20-layer slab model along the (001) direction. While the slab bandstructure is by construction independent of k z , signatures closely related to the bulk k z dispersion are observed in photon energy dependent matrix element calculations, as seen in Fig. 6b. We estimate the attenuation factor e −ξ/2zi of the escaping photoelectrons using the universal escape depth curve from Ref. 27. The k z value probed is calculated using an inner potential of V 0 = 12.2 eV [10]. While at both low and high photon energies the penetration depth is sufficiently large that that ∆k z should be less than 0.05Å −1 , at hν = 71 eV ∆k z = 0.11Å −1 , and linewidth broadening is observed as a result (note here π/c = 0.57Å −1 ). Conversely, for larger ξ values, ∆k z becomes negligible, and something akin to the bulk electronic structure is recovered. Note that ∆k z is not explicitly included in these calculations, but emerges naturally from the combination of slab geometry with variable penetration depth. From these results, it becomes evident that the surface sensitivity can complicate successful estimation of the bulk electronic structure. This emphasizes the need for proper characterization of the k z dispersion, accessible via photon-energy dependent measurements, as in Fig. 6b.

Photoelectron Interference, and Spin-ARPES
Despite these challenges associated with surface sensitivity, the surface can also precipitate new states localized to the interface region which are not possible in bulk systems. Such is the case for example in the Shockley surface states observed along the (111)-termination of noble metals [23,44], Fermi arcs on Weyl-semimetals [45], and conductive surface states observed in topological insulators such as Bi 2 Se 3 [46].
To model the ARPES spectra from these surface states, an extended lattice basis is required, with the unit cell projected onto a slab-geometry. Our implementation of the slab generation is inspired by the algorithm in Ref. 47, allowing for nearly total automation of the slab Hamiltonian initialization. Given a surface Miller index, new lattice vectors can be defined which projects the new unit cell along the desired surface direction to the desired thickness. The bulk Hamiltonian can be propagated over this slab supercell. While we formally maintain periodic boundary conditions, rather than preserve the full translational symmetry of the bulk crystal, a vacuum buffer is defined, with a thickness sufficiently large to suppress hopping elements between neighbouring slab unit cells. The precise location of the crystal-vacuum interface is tuned by the user to achieve the desired surface termination of the crystal. In the case of Bi 2 Se 3 , this termination must occur between the van der Waals-bonded layers of two adjacent quintuple layers (QL) to preserve the topological surface states. The procedure is illustrated in Fig.  7a-c.
Expansion of the basis set to a suitably large slab carries the caveat of a significant memory overhead, which can be to an extent mitigated in the calculation of ARPES intensity: as the finite penetration depth of the probe and photoemitted electrons limit the volume of the unit cell to which we are actually sensitive, the eigenvectors are truncated beyond a modest multiple of the mean free path, allowing for both efficient and highfidelity surface-projected ARPES maps to be computed, as done for the 400-orbital basis used for the simulation in Fig. 8a. As a result, it is the mean-free path more than the size of the basis which limits the ability to treat very large slab unit cells. As an example of this functionality, simulated and experimental ARPES intensity from Bi 2 Se 3 are plotted in Fig. 8. Many of the central tenets of a model strong TI have been confirmed in this material, such as the anticipated chiral spin texture, observed directly via spin-resolved ARPES [46,[48][49][50]. Such spin-resolved experiments [51][52][53] can also be simulated within the chinook software, as shown in Fig. 8b, where we present the simulated spin polarization: This result is in agreement with experiment [46], and can be compared favourably with the surface-projected expectation value of the spinŜ operator, e − |ẑ| ξŜ , plotted in Fig 7e. The bulk states, which lack any discernible spin-polarization (Fig. 7d), vanish from the calculation of P y and so do not appear in Fig. 8b.
While the topological surface states Ψ T SS are primarily composed of p z orbitals at the surface, a pronounced modulation of the photoemission intensity around the Dirac cone is observed as a consequence of the finite extension of Ψ T SS into the crystal bulk. Hybridization with bulk states, in addition to interlayer photoelectron interference can be understood as the progenitor of this modulation, as explored in depth in Ref. 8. The interpretation of this angular intensity pattern in ARPES measurements presents an essential experimental verification and explanation of the limitations of applying a simple k · p model to the description of real topological insulators such as Bi 2 Se 3 . While localized within a finite region near the vacuum interface, the full three-dimensionality of the surface state becomes apparent through consideration of this spectroscopic evidence. Convenient extension to a slab-geometry is then a critical functionality offered by the chinook package.

Variable Experimental Geometry
In practice, ARPES experiments rotate either the sample normal or spectrometer in order to access a broad set of emission angles. While some modern techniques such as photoemission electron microscopy (PEEM) [54], angle-resolved time-of-flight (ARTOF) [48,55], and deflector-based ARPES [56][57][58] apparatus avoid this complication, the assumption of a constant experimental geometry is not always possible. Furthermore, it is often advantageous to rotate the sample orientation in order to for example explore large regions of momentum space, or to achieve better momentum resolution available at higher emission angles [3]. While the most direct complication is associated with variable photon polarization, in the case of S-ARPES, the relative orientation of the detector with respect to the sample is essential to interpret data correctly. To exemplify the practical considerations associated with such experiments wherein the geometry is variable, one may consider an exploration of Rashba-split spin-polarized surface states, as in for example PtCoO 2 [59]. Using the model presented in Ref. 59, in Fig. 9, the Fermi surface is surveyed over several neighbouring Brillouin zones, accessed by rotation of the sample about the horizontal (i.e. k x ) axis. The ARPES intensity calculated with and without consideration of the rotated polarization vector can be compared in Fig. 9a, b. While the intensity in the first Brillouin zone is fairly homogeneous in either case, the higher order zones reflect more substantial variation. These rotations complicate the extraction of orbital character from the photoemission intensity, requiring explicit consideration of the polarization rotation.
In the context of spin-resolved measurements, the spinprojection is measured within the laboratory frame-ofreference, which remains fixed for all sample orientations. As the sample is rotated, contamination of orthogonal In (d-f) we compare the measured spin polarization in the first (d) and second Brillouin zones (e,f). Axis labels indicate distance from Γ1,2. While a chiral Rashba spin texture is observed near normal emission, contamination between the different spin-channels is manifest as a substantial and artificial out-of plane spin projection in the second zone. Pz is zero near Γ1 and not shown here. All related colourscales are represented on the same scale to facilitate direct comparison.
spin-channels is inevitable. While such effects are minimal near Γ 1 , a significant out-of-plane spin-polarization arises near Γ 2 , with P z over 36% of P x , as demonstrated in Fig. 9d-f. It is important to note that the intrinsic spin-polarization is confined entirely to the plane; this apparent out-of-plane polarization exists only in the coordinate frame of the laboratory apparatus. Accounting for the rotation of the measurement coordinate frame for the given experiment, one can then redistribute this information into the channels associated with the sample's intrinsic spin texture: where here P exp is the measured spin polarization from Eq. 18 and S(θ) the spin-projection axis measured at each emission angle. Although Fig. 9b indicates higher photoemission intensity is available in the second Brillouin zone, Fig. 9e,f illustrate the practical challenges associated with resolving the spin-texture near Γ 2 . By affording the user with an ability to encode a realistic experimental configuration in the simulation, such effects can be accounted for in detail, circumventing a signif-icant experimental limitation which may otherwise restrict more general application of the techniques detailed here.

CONCLUSION
We have presented here a simple and powerful numerical framework implemented in Python for the simulation and interpretation of ARPES spectra for a broad variety of materials of interest. Designed with this specific purpose in mind, the open-source structure of the chinook software package is engineered to accommodate further extension beyond this application, as we have done recently for the study of resonant optical excitations in pump-probe spectroscopy experiments [60]. Through the development of these tools, we hope to motivate and facilitate the consideration of the great depth of information encoded in the matrix element of angle-resolved photoemission spectroscopy towards a better understanding of these experiments and the electronic structure of the materials under consideration.  An essential point we want to emphasize is that while a distinction between radial integrals appears to be of critical importance in reproducing experimental data, much qualitative detail seems resilient against the precise relative cross sections. We now make this clear for a simple test case. One of the primary flaws of M F T is its prediction of nodes in the photoemission intensity whenever k ·ˆ = 0. For in-plane polarization then, these effects are anticipated to be quite pronounced. As is seen for the case of FeSe, under M F T , intensity should vanish for k y = 0 in panel (a) of Fig. 1 of the main text. On the contrary, this is where intensity is in fact largest. Taking this as a motivating example, we construct two circular Fermi surface contours defined by helical orbital texture, as anticipated for hole-doped Fe-based superconductors [10,11].
To make these claims more generic, we perform a similar analysis for analogous wavefunctions composed of p, rather than d orbitals. Explicitly, we use the following wavefunctions: Fe: Z = 26 n = 3 l = 2 ψ outer = cos(φ)d yz − sin(φ)d xz ψ inner = sin(φ)d yz + cos(φ)d xz C: Z = 6 n = 2 l = 1 ψ outer = cos(φ)p y − sin(φ)p x ψ inner = sin(φ)p y + cos(φ)p x The subscripts inner (outer) refer to the states with smaller (larger) k F . In addition, φ = tan −1 (k y /k x ). The orbital projections are indicated schematically in the inset of panel (a) of the figure. As elsewhere, we assume a free electron final state in our calculation of the ARPES matrix element.
In the figure, we compare and contrast the effect of identical (B l+1 = B l−1 ) and distinct (B l+1 = B l−1 ) radial integrals on this model system. Forˆ = [0, 1, 0], M F T predicts a node in the photoemission intensity along k y = 0. A reasonable metric then for comparison of B l−1 /B l+1 is the ratio of maximum intensities along the y axis (I y where k x = 0) and x axis (I x , k y = 0), labeled here as I x /I y . The specific momentum points used