# An open-source, end-to-end workflow for multidimensional photoemission
spectroscopy

R. Patrick Xian  Fritz Haber Institute of the Max Planck Society, 14195
Berlin, Germany  Corresponding authors<EMAIL_ADDRESS>rettig@fhi-
berlin.mpg.de<EMAIL_ADDRESS>Yves Acremann  Laboratory for
Solid State Physics, ETH Zurich, 8093 Zurich, Switzerland  Steinn Y. Agustsson
Department of Physics, University of Mainz, 55128 Mainz, Germany  Maciej
Dendzik  Fritz Haber Institute of the Max Planck Society, 14195 Berlin,
Germany  Kevin Bühlmann  Laboratory for Solid State Physics, ETH Zurich, 8093
Zurich, Switzerland  Davide Curcio  Department of Physics and Astronomy,
Interdisciplinary Nanoscience Center (iNANO), Aarhus University, 8000 Aarhus
C, Denmark  Dmytro Kutnyakhov  DESY Photon Science, 22607 Hamburg, Germany
Federico Pressacco  DESY Photon Science, 22607 Hamburg, Germany  Department of
Physics, University of Hamburg, 22761 Hamburg, Germany  Michael Heber  DESY
Photon Science, 22607 Hamburg, Germany  Shuo Dong  Fritz Haber Institute of
the Max Planck Society, 14195 Berlin, Germany  Tommaso Pincelli  Fritz Haber
Institute of the Max Planck Society, 14195 Berlin, Germany  Jure Demsar
Department of Physics, University of Mainz, 55128 Mainz, Germany  Wilfried
Wurth  DESY Photon Science, 22607 Hamburg, Germany  Department of Physics,
University of Hamburg, 22761 Hamburg, Germany  Philip Hofmann  Department of
Physics and Astronomy, Interdisciplinary Nanoscience Center (iNANO), Aarhus
University, 8000 Aarhus C, Denmark  Martin Wolf  Fritz Haber Institute of the
Max Planck Society, 14195 Berlin, Germany  Markus Scheidgen  Fritz Haber
Institute of the Max Planck Society, 14195 Berlin, Germany  Department of
Physics, Humboldt University of Berlin, 12489 Berlin, Germany  Laurenz Rettig
Fritz Haber Institute of the Max Planck Society, 14195 Berlin, Germany
Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Ralph Ernstorfer  Fritz Haber Institute of the
Max Planck Society, 14195 Berlin, Germany  Corresponding authors: xian@fhi-
berlin.mpg.de<EMAIL_ADDRESS><EMAIL_ADDRESS>
> ## Abstract
>
> Characterization of the electronic band structure of solid state materials
> is routinely performed using photoemission spectroscopy. Recent advancements
> in short-wavelength light sources and electron detectors give rise to
> multidimensional photoemission spectroscopy, allowing parallel measurements
> of the electron spectral function simultaneously in energy, two momentum
> components and additional physical parameters with single-event detection
> capability. Efficient processing of the photoelectron event streams at a
> rate of up to tens of megabytes per second will enable rapid band mapping
> for materials characterization. We describe an open-source workflow that
> allows user interaction with billion-count single-electron events in
> photoemission band mapping experiments, compatible with beamlines at
> $3^{\text{rd}}$ and $4^{\text{th}}$ generation light sources and table-top
> laser-based setups. The workflow offers an end-to-end recipe from
> distributed operations on single-event data to structured formats for
> downstream scientific tasks and storage to materials science database
> integration. Both the workflow and processed data can be archived for reuse,
> providing the infrastructure for documenting the provenance and lineage of
> photoemission data for future high-throughput experiments.

## Introduction

Many disciplines in the natural sciences are increasingly dealing with densely
sampled multidimensional datasets. The scientific workflows to obtain and
process them are becoming increasingly complex due to the provenance and
structure of the data and the information needed to be extracted and analyzed
[1, 2]. In materials science and condensed matter physics, various
spectroscopic and structural characterization techniques produce experimental
data of distinct formats and characteristics. Their creation and understanding
require customized processing and analysis pipelines designed by specialists
in the respective fields. The growing incentive for building experimental
materials science databases [3] that complement established theoretical
counterparts [4] calls for open-source and reusable workflows for data
processing [5, 6] that transform raw data to shareable formats for downstream
query, analysis and comparison by non-specialists of the experimental
techniques [7, 8]. Among the various properties associated with materials, the
electronic band structure (EBS) of condensed matter systems is of vital
importance to the understanding of their electronic properties in and out of
equilibrium. Multidimensional photoemission spectroscopy (MPES) [9, 10, 11] is
an emerging technique that bears the potential of high-throughput EBS
characterization through band mapping experiments and holds promise as an
enabling technology for building experimental EBS databases, where data
integration requires traceable knowledge of the processing steps between the
archived and the raw format. Here we present an open-source workflow that
focuses on band mapping data from MPES. In the following, we briefly introduce
the technology of MPES and the associated data processing, before providing
details on the workflow from raw data to database integration.

Figure 1: Schematic of the workflow in MPES. The data acquisition in MPES
starts from (1) photoelectrons liberated by the extreme UV (XUV) or X-ray
photons travelling through the lens systems and the TOF tube to trigger
detection events on the delay-line detector (DLD). (2) Single-event data
acquisition is monitored and controlled by the measurement controller
computer. The raw data are first streamed and stored onto a hard drive in HDF5
format (.h5) and subsequently processed in the workflow through (3) file
reduction (optional), (4)(6) distributed binning, (5) artifact correction and
axis calibrations, carried out at the single-event or the binned data levels.
At the end of the workflow, other data formats are generated (such as HDF5,
MAT or TIFF) for (7) storage, visualization or downstream analysis for
extracting relevant physical parameters. Critical parameters within the
workflow may be exported (as workflow parameters files), shared and reused for
processing other datasets.

MPES, also called momentum microscopy (MM), is born out of the recent
integration of time-of-flight (TOF) electron spectrometers with delay-line
detectors (DLDs) and improved electron-optic lens designs [12, 13, 14, 15].
Compared with the earlier generations of angle-resolved photoemission
spectroscopy (ARPES) [16, 17] using hemispherical analyzers to measure the 2D
energy-momentum distribution of the photoemitted electrons [18], MPES is
capable of recording single-electron events simultaneously sorted into the
$(k_{x},k_{y},E)$ coordinates ($E$: electron energy, $k_{x},k_{y}$: parallel
momentum components) in band mapping experiments, obviating the need for
scanning across sample orientations and subsequent data merging as is the case
for similar experiments using a single hemispherical analyzer. Operation of
the TOF DLD in MPES requires a pulsed photon source and is directly compatible
with 3${}^{\text{rd}}$ and 4${}^{\text{th}}$ generation light sources [19] as
well as laboratory-based table-top setups [20, 21, 22, 23], harnessing their
high repetition rates in the range of multi-kilohertz to megahertz to
drastically improve the detection speed and efficiency. Mapping of the 3D band
structure with sufficient signal-to-noise ratio (SNR) may be achieved on the
timescale of minutes. The technological convergence opens up the possibilities
to record 3D datasets in dependence of one or more additional parameters, such
as spatial location $I(x,y,k_{x},k_{y},E)$, probe photon energy,
$I(k_{x},k_{y},E,k_{z})$ [10], spin-polarization, $I(k_{x},k_{y},E,S)$ [9], or
pump-probe time in time-resolved MM, $I(k_{x},k_{y},E,t)$ [24] within a
reasonable time frame.

From the data perspective, the pulsed sources with high repetition rates
generate densely sampled data at rates of multiple megabytes per second
(MB/s), which has brought about challenges in data processing and management
compared with conventional ARPES experiments. The raw data in MPES are single
photoelectron events registered by the DLD and the physical quantities related
to the detected events are streamed in parallel to the storage files in a
hierarchical file format (e.g. HDF5 [25]). A typical dataset involves
$10^{7}-10^{10}$ detected events with a total size of up to a few hundred of
gigabytes (GBs), depending on the number of coordinates measured (3D or 4D)
and the required SNR. Unlike the large 2D or 3D image-based datasets, such as
those obtained in various forms of optical [26, 27] and electron microscopy
techniques [28, 29], processing and conversion of tabulated single-event data
requires additional steps of statistical computing for conversion into
standard images. This motivates the current workflow development for efficient
data processing and analysis. In data processing and calibration, experiments
performed at different facilities share similar procedures going from the raw
events to the multidimensional hypervolume with calibrated axes, which is the
basis for archiving and downstream analysis. To maintain reproducibility for
the particular data source characteristics, data structure and processing
procedure, we have summarized the workflow (see Fig. 1) into two open-source
software packages (hextof-processor [30] and mpes [31]), with similar design
principles for coping with large-scale facility and table-top experiments,
respectively. The core of our approach includes distributed statistical
processing at the single-event level using parameters calibrated and
determined from preprocessed volumetric datasets, which enables effective
instrument diagnostics, artifact correction, and sample condition monitoring.
The algorithms involved balance physical knowledge and existing methods in
image processing and computer vision. The workflow is illustrated next with
data obtained at some of the electron momentum microscopes currently in
operation, such as the HEXTOF (high energy X-ray time-of-flight) measurement
system [24] at the free-electron laser source FLASH [32] at DESY, and the
table-top high harmonic generation-based setup at the Fritz Haber Institute
(FHI) [21] involving a commercial TOF and DLD (METIS 1000, SPECS GmbH). We use
the material example of tungsten diselenide (WSe2) measured at both
experimental setups to demonstrate the workflow execution, because in momentum
space, the patent features of WSe2 band structure and the nonequilibrium
dynamics initiated by optical excitation of WSe2 have been thoroughly studied
in the past (see Methods) [33, 34, 35, 36, 24]. We expect the workflow
described here to serve as a blueprint for upcoming software platforms in
similar setups to be installed in other facilities or laboratories worldwide.

## Results

  1. Workflow description. The workflow schematic shown in Fig. 1 starts with raw single-event data from measurements. The data are (i) binned in a distributed fashion in the measurement coordinates, including each of the photoelectrons’ position on the detector $(X,Y)$, its TOF, a digital encoder (ENC) axis, and others, if more than four dimensions are acquired in parallel. The binned histogram is (ii) used to estimate the numerical transforms for distortion correction and axis calibration. Next, these transforms are (iii) applied to the raw single-event data to convert the measurement coordinates to the physical axes, $(k_{x},k_{y},E,t_{\text{pp}})$ and others for higher dimensions (see also Fig. 2). Finally, the single-event data are (iv) binned in the transformed grid to yield 3D, 3D+t or other higher-dimensional data with the correct axis values. The outcome may be exported in different formats for storage, visualization and downstream analysis.

  2. Tasks and software infrastructure. Processing billion-count single-event data requires user interaction for data checking and distributed processing to reduce the time consumption. The general tasks in the workflow include the transformation of data streams to multidimensional histograms, artifact correction and axis calibration. These operations can be efficiently decomposed into column-wise operations of the distributed dataframe format offered by the dask package [37] in Python. While the use of dask dataframes provide the common foundation for interactivity with single events of hextof-processor and mpes, they distinguish themselves by the experimental requirements.

Figure 2: Examples of workflow components. Illustrations are given for
artifact correction and axis calibration. Characteristic 1D distributions of
the measured $X$, $Y$, TOF, ENC and an arbitrary axis are shown on the very
left. $U(0,1)$ represent uniformly distributed random noise added to suppress
digitization artifacts (jittering or dithering). The transforms ($g$’s) are
calibration functions that convert the values in the measurement axes to the
physical ones. The transform $\mathcal{L}(X,Y)$ corrects the symmetry
distortion, while the spherical timing aberration and space charge are
compensated for by $\Delta\text{TOF}_{\text{sph}}$ and
$\Delta\text{TOF}_{\text{sc}}$, respectively. Binning of the corrected single-
event data over the calibrated physical axes yields a multidimensional
hypervolume (right picture) of photoemission intensity data along with the
physical axes values.

At large-scale facilities, experiments often record a large number of machine
parameters that need to be stored, though only a small number of relevant
parameters are needed for downstream processing. Therefore, the hextof-
processor package includes a parameter sampling step to retrieve intermediate
tabulated data in the Apache Parquet format (https://parquet.apache.org/), a
column-based data structure optimized for computational efficiency. This
approach reduces the processing overhead in searching through the raw data
files every time when data are queried during the subsequent processing. As an
open-source project, other beamtime-specific functionalities are added by
users to the existing framework at every new experimental run. The mpes
package adapts to the much simpler file structure produced at table-top
experimental setups and makes direct use of the HDF5 raw data. It comes with
added functionalities motivated by the existing issues encountered in data
acquisition and downstream processing such as axis calibration, masking,
alignment and different forms of artifact correction. The softwares come with
detailed documentation and examples online for users to gain familiarity (see
Code Availability).

  3. Artifact correction. Artifacts in MPES data come from mechanical imperfections, stray fields (electric and magnetic), uncertainties in the alignment of the sample, light beams and the multistage electron-optic lens systems as well as the data digitization process. Minimizing and correcting instrumental imperfections plays an important role in the validity of downstream analysis. We carry out artifact correction sequentially at the level of single photoelectron events or the data hypervolume obtained from multidimensional histogramming (see Fig. 2). The outcomes are illustrated using the correction of (1) digitization artifact (see Fig. 3) and (2) spherical timing aberration artifact (see Fig. 4), with technical details in Methods.

Figure 3: Digitization artifact correction by histogram jittering. Removal of
the digitization artifact is illustrated with a 2D $k$-$E$ slice across the
Brillouin zone center (at 0 $\mathrm{\AA}^{-1}$ momentum) of the band mapping
dataset measured at FHI on WSe2. The images before and after histogram
jittering and their difference are shown in a, b and c respectively. A zoomed-
in section of the data are shown in the insets in a-c. The effective removal
of the digitization artifact is further demonstrated in the momentum-
integrated energy distribution curves in d. The traces in d are computed by
averaging horizontally over their corresponding 2D images in a-c.

  4. Axis calibration. To transform the measurement axes of the DLD into physically relevant axes for electronic band mapping, calibrations are required, as shown in Fig. 2. The calibration functions are constructed with parameters derived from comparing physical knowledge of the materials (e.g. Brillouin zone size, Fermi level position) with the corresponding scales in data. They are applied either to the binned data hypervolume, or to the single-electron events raw data individually in a distributed fashion before binning. Details on the calibration data transforms are provided in Methods.

Figure 4: Spherical timing aberration correction. The correction is
demonstrated using W4f core-level data measured at FLASH. The energy spacing
between the W4f7/2 and W4f5/2 levels is about 2.1 eV [34]. a. Illustration of
the geometric origin of the spherical timing aberration in the time-of-flight
(TOF) drift tube. b. Comparison of the W4f spectra at the center and on the
edge of the detector plane. The energy spectra are extracted from the
corresponding regions, marked by the dots in the same blue and red colors,
respectively, in c. The white stripes crossing at the detector center block
the exposed edges of the four-quadrant detector quadrants. d. The uncorrected
and corrected radial-averaged peak TOF positions for the W4f7/2 core level.

  5. Data storage and format. The simplistic form of the output data hypervolume derived from single-electron events includes non-negative scalar values of the photoemission intensity and the calibrated real-valued axes coordinates, including $k_{x}$, $k_{y}$, $E$, and other parameter dependencies such as the pump-probe time delay $t_{pp}$. These values are exported as HDF5, MAT or TIFF, with the metadata included as attributes of the files.

  6. Workflow archiving and reuse. Computational workflows are valued by their the reproducibility [38]. Archiving and sharing the workflow parameters among users of the beamlines or facilities allow comparison between experimental runs and reuse for the simultaneous benefits of machine diagnostics and experimental efficiency. To achieve this, we store critical parameters generated within the workflow in a separate file as workflow parameters (see Fig. 1) during each step, including the numerical values used in binning, the intermediate parameters and coefficients of the correction and calibration functions, etc. They can be reused when loading into the processing of other datasets.

Figure 5: Typical visual representations of the volumetric band mapping data.
The examples are illustrated using band mapping data of the layered
semiconductor WSe2, measured with the HEXTOF instrument at FLASH (a,c) and the
METIS detector at the FHI (b). The visualizations are a. the orthoslice
representation, b. the band-path diagram (right) with the momentum path
labelled in dashed blue line in the momentum $k_{x}$-$k_{y}$ plane (left), and
c. the cut-out view. All color scales represent photoemission intensity. The
letters label the high symmetry points of the hexagonal Brillouin zone of WSe2
[33]. The intensity features on the upper side nearby the
$\overline{\mathrm{K}}$ and $\overline{\mathrm{K}^{\prime}}$ points result
from nonequilibrium electronic populations following photoexcitation of
electrons to the conduction band, while the intensities below are from valence
band photoemission (see Methods for experimental conditions).

  7. Data visualization. The adaptation of established scientific visualization methods in the physical sciences [39, 40] to band mapping data should incorporate the requirements and knowledge of the data characteristics in this field of research. The band mapping data in 3D (multi-megavoxel) and 3D+t (multi-gigavoxel) include the inherent symmetries from the electronic band structure of the material, but the intensity modulations in the photoemission process [41], dynamics and sample condition disrupt the original symmetry. The overall goal is to emphasize the features of interest while exploiting the symmetry to simplify the visualization (see Methods). The output files from the processing pipeline are compatible with open-source visualization software such as matplotlib [42], ParaView [39] and Blender [43].

  8. Downstream analysis integration. Typical photoemission data analysis involves extracting electronic band structure parameters, physical coupling constants and lifetimes via fitting of lineshapes [16] or dynamical models [44], often carried out specific to the material under study. At the end of our distributed workflow, the data size is on the order of a few to tens of gigabytes, which can be directly loaded into memory on users’ local machines for downstream data analysis with custom routines.

  9. Experimental metadata. The metadata of the data files have a tree structure and contain information of the experimental setting, parameters of the pulsed light source, the detector and the sample under study. A list of top-level metadata parameters is presented in Table 1. A full and current list of all metadata parameters, including the top-level parameters and their constituent lower-level parameters, along with their definitions, units and values, is provided in Supplementary Tables 1-4. For database integration, an accompanying data parser (parser-mpes, see Code Availability) for MPES data has been written in accordance with existing standards [45] for computational materials science in NOMAD [8], featuring an electronic version of the metadata parameter list in the file mpes.nomadmetainfo.json online. The metadata parameter list and the data parser are versioned and are updated based on the corresponding changes in the data structure for photoemission spectroscopy experiments. The existing WSe2 photoemission data have been integrated into the experimental section of the materials science database NOMAD (see Data Availability).

Table 1: Top-level metadata parameters Category name | Description
---|---
General parameters | Descriptive information of the experiment and facility
Source parameters | Technical parameters relating to the photon source
Detector parameters | Technical parameters relating to the photoelectron detector
Sample parameters | Parameters relating to the material sample in experiment

## Discussion

We have designed and implemented an open-source, end-to-end workflow for
processing single-event data produced in multidimensional photoemission
spectroscopy, linking to downstream tasks, providing guidelines and software
for integrating processed data into the NOMAD experimental materials science
database. The distributed processing takes full advantage of the single-event
data streams directly accessible from the TOF delay-line detector for event-
wise correction and calibration and converts the raw events to the calibrated
data hypervolume for project-specific downstream analysis. The functionalities
within the workflow are publicly accessible through the software packages we
have developed (hextof-processor [30] and mpes [31]). The processing workflow
is archived at each step of operation and the processed data may be integrated
into experimental database with user-specified metadata. The methods described
here are applicable to all existing types of multidimensional photoemission
band mapping measurements beyond the static and time-dependent settings
described here.

Our end-to-end workflow from raw data to processed data to database
integration provides a fast-track and all-in-one solution to the demands for
open experimental data and reproducible research in the materials science
community [7, 8]. The public repositories for the software packages are the
foundations for phased future extension and integration with existing
analytical tools in the photoemission spectroscopy community. The modular
structure of the packages introduced here allows targeted upgrades by both
temporary and dedicated maintainers and users. Casting the workflow in the
Python programming environment provides the foundation for convenient
incorporation of existing image processing and machine-learning resources [46]
for further exploration and understanding of the band mapping datasets, which
contain rich information owing to the complex nature of the photoemission
process [16, 18]. This is especially beneficial for broader adoption of
photoemission since the interpretation of photoemission data is often linked
to the observed or extracted outstanding features such as local intensity
extrema, dispersion kinks and satellites, lineshape parameters and pattern
symmetry [16], therefore, the access to experimental data and the potential
integration with existing electronic structure-related software [5, 47, 48,
49] will facilitate method developments and the direct comparison between
experimental results and theoretical band structure calculations within the
same programming platform.

## Methods

  1. Sample preparation. Single-crystalline samples of 2D bulk WSe2 (2H stacking) were purchased from HQ Graphene. Crystals of size around 5 mm $\times$ 5 mm $\times$ 1 mm were used directly for the measurements. To prepare a clean surface by cleaving, we attached a cleaving pin upright to the sample surface using conducting epoxy (EPOTEC H20) outside the vacuum chamber and removed the pin by mechanical force in ultrahigh vacuum.

  2. Photoemission experiments. The measurements were conducted using the HEXTOF instrument [24] at the DESY FLASH PG-2 beamline [50] with the free-electron laser (FEL) as well as a laboratory source [21] with a METIS electron momentum microscope (SPECS METIS 1000) installed at the FHI. In the measurements at FLASH, the FEL was tuned to 36.5 eV (or 34.0 nm) and 109 eV (or 11.4 nm), the optical pump pulse had a center wavelength of 775 nm. The measurements at the FHI used a 21.7 eV home-built extreme UV source based on high harmonic generation in Ar gas driven by an optical parametric chirped-pulse amplifier operating at 500 kHz repetition rate [51]. The optical pump pulse is centered at 800 nm. In both FEL and laboratory experiments, the near-infrared light pulses promote the electronic population at the K and K′ high-symmetry points (corresponding to $\overline{\mathrm{K}}$ and $\overline{\mathrm{K}^{\prime}}$ points, respectively, in the projected Brillouin zone obtained from photoemission, as shown in Fig. 5) in momentum space to the excited states via direct optical transitions. The nonequilibrium electronic dynamics are probed via valence and conduction band photoemission [35] as well as core-level photoemission [36], using $s$-polarized extreme UV and soft X-ray probe pulses, respectively.

  3. Digitization artifact. The time-to-digital converter (TDC) outputs digitized data according to the binning width of the on-board electronics. Data conversion from one digitized format to another in a rebinning process often creates a picket fence-like effect (see Fig. 3). This phenomenon originates from the incommensurate bin size in the two rounds of sampling processes (binning and rebinning). To solve the problem, one introduces a slight amount of uniformly distributed noise, with an amplitude equal to half of the original bin size, to the single-event values when carrying out the bin counts. This is similar to the histogram jittering (or dithering) technique [52, 53] used in statistical visualization and computer graphics. Mathematically, the uniformly distributed noise $U(0,1)$ bounded in the range $[0,1]$ is added before binning a univariate data stream, $S=\\{S_{i}\\}$ via,

$S_{i}^{\prime}=S_{i}+\frac{w_{b}}{2}\times U(0,1).$ (1)

Here, $w_{b}$ is the bin width. For binning of multivariate data streams, such
as the detector X position (or $k_{x}$), Y position (or $k_{y}$), and the
photoelectron TOF (or $E$), we adopt the same approach individually for each
dimension. The effect of jittering in reducing the digitization artifact is
demonstrated in Fig. 3.

  4. Spherical timing aberration. Electrons entering the TOF tube at different lateral positions travel through different path lengths to reach the detector, which is the origin of the spherical timing aberration as illustrated in Fig. 4. The lateral position-dependent time delay may be expressed as,

$\Delta\text{TOF}_{\text{sph}}(r)=\big{(}\sqrt{1+r^{2}/d^{2}}-1\big{)}\text{TOF}_{0},$
(2)

where $r$ is the radial distance from the center of the DLD and
$\text{TOF}_{0}$ is the TOF normalization constant. For a typical field-free
region length of $d\sim 1$ m in the TOF tube and a DLD screen radius of $r=50$
mm, $\Delta\text{TOF}/\text{TOF}_{0}\approx 1.25\times 10^{-3}$. Assuming
$\text{TOF}_{0}=0.5$ $\mu s$, the spherical timing aberration in TOF scale is
$\Delta\text{TOF}_{\text{sph}}\approx 0.62$ ns, which is larger than the DLD’s
temporal resolution of $\sim$ 0.15 ns. The effect of the spherical timing
aberration is visible for a few eV energy range with fine bins but quite small
on a large energy range. To illustrate this effect, we use the W$4f$ core-
level data presented in Fig. 4b. For every $(X,Y)$ position on the detector
the peak of W$4f_{7/2}$ was fitted with a Voigt profile and the peak positions
are shown in Fig. 4c. As the spectra from deep core levels typically do not
show dispersion, the deviation from fitting corresponds to the spherical
timing aberration of the electron optics. In order to compensate for the
spherical timing aberration, we first transform the data from Cartesian to the
polar coordinates (see Fig. 4c), and then fit the radial-averaged peak
position to a polynomial function of the radius,

$\Delta\text{TOF}_{\text{sph}}(r)=\frac{r^{2}\text{TOF}_{0}}{2d^{2}}-\frac{r^{4}\text{TOF}_{0}}{8d^{4}}+\text{O}(r^{6}).$
(3)

The fitting results together with the corrected radial distribution are
presented in Fig. 4d.

  5. Symmetry distortion. Photoemission patterns in the $(k_{x},k_{y})$ plane (i.e. an energy slice) may exhibit distorted symmetry due to the influence of various factors from the instrument, the sample and the experimental geometry on the trajectory of low-energy photoelectrons. Correction of the symmetry distortion yet preserving the intensity features requires the use of symmetry-related landmarks to solve for the symmetrization coordinate transform in the framework of nonrigid image registration [54]. In typical situations with an excellent electron lens alignment, the energy dependence of the momentum distortion within the focused phase space volume covering an energy range of several eV is negligible, so the same coordinate transform can be applied to all energy slices in the volumetric data (including both valence and conduction bands) or simultaneously to all single events.

  6. Other single-experiment artifacts. (1) Momentum center shift: The momentum center of the emergent photoelectrons travelling through the electron-optic system may experience an energy-dependent shift owing to the slight misalignment in the system or the influence of stray fields. Correction of the center shift requires an energy-dependent center alignment of energy slices. The shift along the energy (or TOF) axis may be estimated using phase correlation [55] or mutual information-based [56] sequential image registration methods, in which the series of energy slices are treated as an image sequence. In a well-shielded and well-aligned electron-optic lens system, generally, the momentum center shift is negligible in the focused photoelectron energy range. (2) Space-charge effect (SCE): The secondary photoelectron clouds originating from the probe and pump pulses cause a “doming effect” of the photoemission intensity distribution around the momentum center of the band structure. This is especially visible in systems with a clear Fermi edge [9, 11] or non-dispersing shallow core levels, which may be used as references for calibrating the parameters used for the flattening transform by applying a momentum-dependent shift $\Delta\text{TOF}_{\text{sc}}(k_{x},k_{y})$ in the TOF (or the calibrated energy) coordinate of the single-event data.

  7. Momentum calibration. The scaling factors for momentum calibration are computed by comparing the positions of known high symmetry points in the band structure with their corresponding locations in an energy slice. Suppose $A$ and $B$ are two high symmetry points identifiable (e.g. as local extrema) from the experimental data with pixel positions ($X_{A}$, $Y_{A}$) and ($X_{B}$, $Y_{B}$), and momentum positions, ($k_{x}^{A}$, $k_{y}^{A}$) and ($k_{x}^{B}$, $k_{y}^{B}$), respectively. We calculate the pixel-to-momentum scaling ratios, $f_{X}$ and $f_{Y}$, along the $X$ (column) and $Y$ (row) directions of a 2D $k$-space image, respectively. Then, the momentum coordinate ($k_{x}$, $k_{y}$) at each pixel position ($X$, $Y$) may be derived.

$\displaystyle f_{D}$ $\displaystyle=(k_{d}^{A}-k_{d}^{B})/(D_{A}-D_{B})$ (4)
$\displaystyle k_{d}$ $\displaystyle=f_{D}\times(D-D_{A})\quad(D,d=X,x\text{
or }Y,y)$ (5)

  8. Energy calibration. The calibration requires a set of band mapping data measured at different bias voltages (applied between the material sample and the ground), usually sampled with a spacing of 0.5 V in a range of $\pm$3-5 V around the normally applied bias voltage for a particular sample. The calibration proceeds by finding the TOF feature (e.g. local extrema) correspondences in the 1D energy distribution curves (EDCs) at different biases using the dynamic time warping algorithm [57]. The transformation from the TOF to the photoelectron energy $E$ is approximated as a polynomial function,

$E(\text{TOF})=\sum_{i=0}^{n}a_{i}\text{TOF}^{i}$ (6)

The approximation is sufficiently accurate within a range of $\sim$ 20 eV,
sufficient to cover the entire valence band and some low-lying conduction
bands of typical materials. The polynomial coefficients are determined using
nonlinear least squares by solving $\Delta T\cdot\textbf{a}=\Delta E$, in
which $\textbf{a}=(a_{1},a_{2},...)^{T}$ is the coefficient vector while the
constant offset $a_{0}$ is determined by manual alignment to an energy
reference, such as the Fermi level or valence band maximum. The vector $\Delta
E$ and the matrix $\Delta T$ contain, respectively, the pairwise differences
of the bias voltages and the polynomial terms of differential TOF values. To
calibrate a large energy range including multiple core levels, a piecewise
polynomial may be used [11].

  9. Pump-probe delay calibration. The time origin (“time zero”) in time-resolved photoemission spectroscopy, i.e. the temporal overlap of pump and probe pulses, is determined by fitting of a characteristic trace extracted from the data. Since the readings of the digital encoder (see Fig. 2) are sampled linearly, equally-spaced pump-probe delays are directly convertible from the readings using linear interpolation, given the boundary values of the translation stage positions and the corresponding delay times. For unequally-spaced delays, a delay marker is first added to each data point as a separate column after data acquisition to group together the encoder reading ranges that correspond to the specific time delays. The data binning is carried out over the delay marker column instead of the equally-sampled encoder readings.

  10. Visualization strategies. We discuss here three methods for the display of volumetric band mapping data, which are, at the same time, the basis for visualizing 3D+t data with time as an animated axis. (1) The orthoslice representation includes orthogonal 2D planes selected in specific regions in the 3D volume [39], which highlights specific slices deep within the data less visible in a volumetrically rendered view (see Fig. 5a). Along this line, we have developed a software package, 4Dview [58], to explore 4D data using simultaneously linked orthoslices, which also features contrast adjustment and data integration within a hypervolume of interest. (2) The band-path plot (see Fig. 5b) is a 2D representation of the 3D band mapping volume generated by combining a series of 2D cuts along selected momentum paths (or k-paths) traversing a list of so-called high-symmetry points [59, 60]. This representation captures the largest dispersion within the band structure. For volumetric data, the same path may be sampled from all the full energy range to produce the plot shown in Fig. 5b. The analysis and visualization modules in the mpes package include functionalities to compose customized band-path plots. (3) The cut-out view (see Fig. 5c) exposes a specific part of interest in the volumetric data, while not losing the rest [39]. The analysis module in the mpes package provides ways to generate precise cut-outs using position landmarks (e.g. high-symmetry points labelled in Fig. 5) and inequalities.

## Acknowledgements

We thank G. Schönhense for support on the photoelectron detector, S.
Grunewald, S. Schülke and G. Schnapka for support on the computing
infrastructures. We thank G. Brenner, H. Redlin and S. Dziarzhytski at FLASH,
DESY, and H. Meyer and S. Gieschen from the University of Hamburg for beamline
and instrumentation support. The work was partially supported by BiGmax, the
Max Planck Society’s Research Network on Big-Data-Driven Materials-Science,
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (Grant No. ERC-2015-CoG-682843), and the
German Research Foundation (DFG) through the Emmy Noether program under grant
number RE 3977/1 and the SFB/TRR 227 ”Ultrafast Spin Dynamics” (projects A09
and B07). D. Kutnyakhov, M. Heber and W. Wurth acknowledge funding by the DFG
within the framework of the Collaborative Research Centre SFB 925 - 170620586
(project B2). F. Pressacco acknowledges funding from the excellence cluster
EXC 1074 “The Hamburg Centre for Ultrafast Imaging - Structure, Dynamics and
Control of Matter at the Atomic Scale” of the DFG. S. Y. Agustsson and J.
Demsar acknowledge the financial support by the DFG in the framework of the
Collaborative Research Centre SFB TRR 173 - 268565370 (project A5). D. Curcio
and P. Hofmann acknowledge funding from VILLUM FONDEN via the Centre of
Excellence for Dirac Materials (Grant No. 11744). T. Pincelli thanks the
Alexander von Humboldt Foundation for financial support.

## Author contributions

Y.A., K.B., S.Y.A., D.C., R.P.X. and M.D. wrote the hextof-processor package.
R.P.X. and L.R. wrote the mpes package. D.K., Y.A., F.P., R.P.X., S.Y.A.,
D.C., M.D., M.H., S.D., P.H., L.R., R.E. and W.W. participated in the
experiments at the FLASH PG-2 beamline using the HEXTOF instrument in Hamburg.
S.D. and L.R. conducted the experiment at the Fritz Haber Institute using the
METIS electron momentum microscope. R.P.X., M.D., L.R., R.E., M.S., T.P.
constructed the metadata format, R.P.X. and M.S. implemented them into parser-
mpes. R.P.X. wrote the initial manuscript with contributions from M.D. and
Y.A. All authors contributed to discussions to bring the manuscript to its
final form.

## Data availability

The single-event photoemission data used for demonstrating the workflow is
available on the Zenodo platform at 10.5281/zenodo.2704787 (valence and
conduction band photoemission at FEL) [61], 10.5281/zenodo.3945432 (core-level
photoemission at FEL) [62] and 10.5281/zenodo.3987303 (valence band
photoemission from laboratory setup). [63] The preprocessed data are being
integrated into the NOMAD database in the domain for experimental materials
science data accessible at https://nomad-
lab.eu/prod/rae/gui/search?domain=ems.

## Code availability

The code, including documentation and examples in Jupyter notebooks for
implementing the data transformations in the workflow are available as hextof-
processor
(https://github.com/momentoscope/hextof-processor) [30] and mpes
(https://github.com/mpes-kit/mpes) [31]. The parser for integrating
preprocessed experimental data into the NOMAD database is available as parser-
mpes (https://gitlab.mpcdf.mpg.de/rpx/parser-mpes) [64].

## Competing interests

The authors declare no competing interests.

## References

  * [1] Pruneau, C.  _Data Analysis Techniques for Physical Scientists_ (Cambridge University Press, 2017).
  * [2] Deelman, E. _et al._ The future of scientific workflows.  _The International Journal of High Performance Computing Applications_ 32, 159–175 (2018).
  * [3] Zakutayev, A. _et al._ An open experimental database for exploring inorganic materials.  _Scientific Data_ 5, 180053 (2018).
  * [4] Himanen, L., Geurts, A., Foster, A. S. & Rinke, P.  Data‐Driven Materials Science: Status, Challenges, and Perspectives.  _Advanced Science_ 1900808 (2019).
  * [5] Pizzi, G., Togo, A. & Kozinsky, B.  Provenance, workflows, and crystallographic tools in materials science: AiiDA, spglib, and seekpath.  _MRS Bulletin_ 43, 696–702 (2018).
  * [6] Perkel, J. M.  Workflow systems turn raw data into scientific knowledge.  _Nature_ 573, 149–150 (2019).
  * [7] Hill, J. _et al._ Materials science with large-scale data and informatics: Unlocking new opportunities.  _MRS Bulletin_ 41, 399–409 (2016).
  * [8] Draxl, C. & Scheffler, M.  NOMAD: The FAIR concept for big data-driven materials science.  _MRS Bulletin_ 43, 676–682 (2018).
  * [9] Schönhense, G., Medjanik, K. & Elmers, H.-J.  Space-, time- and spin-resolved photoemission.  _Journal of Electron Spectroscopy and Related Phenomena_ 200, 94–118 (2015).
  * [10] Medjanik, K. _et al._ Direct 3D mapping of the Fermi surface and Fermi velocity.  _Nature Materials_ 16, 615–621 (2017).
  * [11] Schönhense, B. _et al._ Multidimensional photoemission spectroscopy—the space-charge limit.  _New Journal of Physics_ 20, 033004 (2018).
  * [12] Krömker, B. _et al._ Development of a momentum microscope for time resolved band structure imaging.  _Review of Scientific Instruments_ 79, 053702 (2008).
  * [13] Ovsyannikov, R. _et al._ Principles and operation of a new type of electron spectrometer – ArTOF.  _Journal of Electron Spectroscopy and Related Phenomena_ 191, 92–103 (2013).
  * [14] Damm, A. _et al._ Application of a time-of-flight spectrometer with delay-line detector for time- and angle-resolved two-photon photoemission.  _Journal of Electron Spectroscopy and Related Phenomena_ 202, 74–80 (2015).
  * [15] Tusche, C., Krasyuk, A. & Kirschner, J.  Spin resolved bandstructure imaging with a high resolution momentum microscope.  _Ultramicroscopy_ 159, 520–529 (2015).
  * [16] Damascelli, A., Hussain, Z. & Shen, Z.-X.  Angle-resolved photoemission studies of the cuprate superconductors.  _Reviews of Modern Physics_ 75, 473–541 (2003).
  * [17] Yang, H. _et al._ Visualizing electronic structures of quantum materials by angle-resolved photoemission spectroscopy.  _Nature Reviews Materials_ 3, 341–353 (2018).
  * [18] Suga, S. & Sekiyama, A.  _Photoelectron Spectroscopy: Bulk and Surface Electronic Structures_ (Springer, 2014).
  * [19] Couprie, M.  New generation of light sources: Present and future.  _Journal of Electron Spectroscopy and Related Phenomena_ 196, 3–13 (2014).
  * [20] Chiang, C.-T. _et al._ Boosting laboratory photoelectron spectroscopy by megahertz high-order harmonics.  _New Journal of Physics_ 17, 013035 (2015).
  * [21] Puppin, M. _et al._ Time- and angle-resolved photoemission spectroscopy of solids in the extreme ultraviolet at 500 kHz repetition rate.  _Review of Scientific Instruments_ 90, 023104 (2019).
  * [22] Corder, C. _et al._ Ultrafast extreme ultraviolet photoemission without space charge.  _Structural Dynamics_ 5, 054301 (2018).
  * [23] Buss, J. H. _et al._ A setup for extreme-ultraviolet ultrafast angle-resolved photoelectron spectroscopy at 50-kHz repetition rate.  _Review of Scientific Instruments_ 90, 023105 (2019).
  * [24] Kutnyakhov, D. _et al._ Time- and momentum-resolved photoemission studies using time-of-flight momentum microscopy at a free-electron laser.  _Review of Scientific Instruments_ 91, 013109 (2020).
  * [25] Folk, M., Heber, G., Koziol, Q., Pourmal, E. & Robinson, D.  An overview of the HDF5 technology suite and its applications.  In _Proceedings of the EDBT/ICDT 2011 Workshop on Array Databases - AD ’11_ , 36–47 (ACM Press, New York, New York, USA, 2011).
  * [26] Weiler, N. C., Collman, F., Vogelstein, J. T., Burns, R. & Smith, S. J.  Synaptic molecular imaging in spared and deprived columns of mouse barrel cortex with array tomography.  _Scientific Data_ 1, 140046 (2014).
  * [27] Ker, D. F. E. _et al._ Phase contrast time-lapse microscopy datasets with automated and manual cell tracking annotations.  _Scientific Data_ 5, 180237 (2018).
  * [28] Levin, B. D. _et al._ Nanomaterial datasets to advance tomography in scanning transmission electron microscopy.  _Scientific Data_ 3, 160041 (2016).
  * [29] Aversa, R., Modarres, M. H., Cozzini, S., Ciancio, R. & Chiusole, A.  The first annotated set of scanning electron microscopy images for nanoscience.  _Scientific Data_ 5, 180172 (2018).
  * [30] Acremann, Y. _et al._ hextof-processor.  https://github.com/momentoscope/hextof-processor (2020).
  * [31] Xian, R. P. & Rettig, L.  mpes.  https://github.com/mpes-kit/mpes (2020).
  * [32] Ackermann, W. _et al._ Operation of a free-electron laser from the extreme ultraviolet to the water window.  _Nature Photonics_ 1, 336–342 (2007).
  * [33] Riley, J. M. _et al._ Direct observation of spin-polarized bulk bands in an inversion-symmetric semiconductor.  _Nature Physics_ 10, 835–839 (2014).
  * [34] Shallenberger, J. R.  2D tungsten diselenide analyzed by XPS.  _Surface Science Spectra_ 25, 014001 (2018).
  * [35] Bertoni, R. _et al._ Generation and Evolution of Spin-, Valley-, and Layer-Polarized Excited Carriers in Inversion-Symmetric WSe2.  _Physical Review Letters_ 117, 277201 (2016).
  * [36] Dendzik, M. _et al._ Observation of an Excitonic Mott Transition Through Ultrafast Core- cum -Conduction Photoemission Spectroscopy.  _Physical Review Letters_ 125, 096401 (2020).
  * [37] Dask Development Team.  _Dask: Library for dynamic task scheduling_ https://dask.org (2016).
  * [38] Stodden, V. _et al._ Enhancing reproducibility for computational methods.  _Science_ 354, 1240–1241 (2016).
  * [39] Hansen, C. D. & Johnson, C. R. (eds.) _The Visualization Handbook_ (Elsevier Butterworth-Heinemann, 2005).
  * [40] Lipşa, D. R. _et al._ Visualization for the Physical Sciences.  _Computer Graphics Forum_ 31, 2317–2347 (2012).
  * [41] Moser, S.  An experimentalist’s guide to the matrix element in angle resolved photoemission.  _Journal of Electron Spectroscopy and Related Phenomena_ 214, 29–52 (2017).
  * [42] Hunter, J. D.  Matplotlib: A 2d graphics environment.  _Computing in Science & Engineering_ 9, 90–95 (2007).
  * [43] Community, B. O.  _Blender - a 3D modelling and rendering package_.  Blender Foundation, Stichting Blender Foundation, Amsterdam  http://www.blender.org (2018).
  * [44] Weinelt, M.  Time-resolved two-photon photoemission from metal surfaces.  _Journal of Physics: Condensed Matter_ 14, R1099–R1141 (2002).
  * [45] Ghiringhelli, L. M. _et al._ Towards efficient data exchange and sharing for big-data driven materials science: metadata and data formats.  _npj Computational Materials_ 3, 46 (2017).
  * [46] Butler, K. T., Davies, D. W., Cartwright, H., Isayev, O. & Walsh, A.  Machine learning for molecular and materials science.  _Nature_ 559, 547–555 (2018).
  * [47] Ong, S. P. _et al._ Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis.  _Computational Materials Science_ 68, 314–319 (2013).
  * [48] Hjorth Larsen, A. _et al._ The atomic simulation environment—a Python library for working with atoms.  _Journal of Physics: Condensed Matter_ 29, 273002 (2017).
  * [49] M. Ganose, A., J. Jackson, A. & O. Scanlon, D.  sumo: Command-line tools for plotting and analysis of periodic ab initio calculations.  _Journal of Open Source Software_ 3, 717 (2018).
  * [50] Gerasimova, N., Dziarzhytski, S. & Feldhaus, J.  The monochromator beamline at FLASH: performance, capabilities and upgrade plans.  _Journal of Modern Optics_ 58, 1480–1485 (2011).
  * [51] Puppin, M. _et al._ 500 kHz OPCPA delivering tunable sub-20 fs pulses with 15 W average power based on an all-ytterbium laser.  _Optics Express_ 23, 1491 (2015).
  * [52] Chambers, M., Cleveland, S., Tukey, A. & Kleiner, B.  _Graphical Methods for Data Analysis_ (Wadsworth International Group, 1983).
  * [53] Novo, D. & Wood, J.  Flow cytometry histograms: Transformations, resolution, and display.  _Cytometry Part A_ 73A, 685–692 (2008).
  * [54] Xian, R. P., Rettig, L. & Ernstorfer, R.  Symmetry-guided nonrigid registration: The case for distortion correction in multidimensional photoemission spectroscopy.  _Ultramicroscopy_ 202, 133–139 (2019).
  * [55] Guizar-Sicairos, M., Thurman, S. T. & Fienup, J. R.  Efficient subpixel image registration algorithms.  _Optics Letters_ 33, 156 (2008).
  * [56] P. Viola & W.M.Wells.  Alignment by Maximisation of Mutual Information.  _International Journal of Computer Vision_ 24, 137–154 (1997).
  * [57] Salvador, S. & Chan, P.  Toward accurate dynamic time warping in linear time and space.  _Intelligent Data Analysis_ 11, 561–580 (2007).
  * [58] Dendzik, M.  4Dview 1.0.  _Zenodo_ https://doi.org/10.5281/zenodo.3360817 (2019).
  * [59] Setyawan, W. & Curtarolo, S.  High-throughput electronic band structure calculations: Challenges and tools.  _Computational Materials Science_ 49, 299–312 (2010).
  * [60] Hinuma, Y., Pizzi, G., Kumagai, Y., Oba, F. & Tanaka, I.  Band structure diagram paths based on crystallography.  _Computational Materials Science_ 128, 140–184 (2017).
  * [61] Xian, R. P. _et al._ Multidimensional photoemission spectra of tungsten diselenide.  _Zenodo_ https://doi.org/10.5281/zenodo.2704787 (2020).
  * [62] Dendzik, M. _et al._ Time-resolved core-level photoemission data of tungsten diselenide.  _Zenodo_ https://doi.org/10.5281/zenodo.3945432 (2020).
  * [63] Xian, R. P. _et al._ Datasets for the computational workflow of multidimensional photoemission spectroscopy.  _Zenodo_ https://doi.org/10.5281/zenodo.3987303 (2020).
  * [64] Xian, R. P. & Scheidgen, M.  parser-mpes.  https://gitlab.mpcdf.mpg.de/rpx/parser-mpes (2019).


# Applicability and limitations of cluster perturbation theory for Hubbard
models

Nicklas Enenkel  HQS Quantum Simulations GmbH, Rintheimer Straße 23, 76131
Karlsruhe, Germany  Institut für Theoretische Festkörperphysik, Karlsruhe
Institute of Technology, 76131 Karlsruhe, Germany  Markus Garst  Institut für
Theoretische Festkörperphysik, Karlsruhe Institute of Technology, 76131
Karlsruhe, Germany  Institute for Quantum Materials and Technology, Karlsruhe
Institute of Technology, 76131 Karlsruhe, Germany  Peter Schmitteckert  HQS
Quantum Simulations GmbH, Rintheimer Straße 23, 76131 Karlsruhe, Germany

###### Abstract

We present important use cases and limitations when considering results
obtained from Cluster Perturbation Theory (CPT). CPT combines the solutions of
small individual clusters of an infinite lattice system with the Bloch theory
of conventional band theory in order to provide an approximation for the
Green’s function in the thermodynamic limit. To this end we are investigating
single-band and multi-band Hubbard models in one- and two-dimensional systems.
A special interest is taken in the supposed pseudo gap regime of the two-
dimensional square lattice at half filling and intermediate interaction
strength ($U\leq 3t$) as well as the metal-insulator transition. We point out
that the finite-size level spacing of the cluster limits the resolution of
spectral features within CPT. This restricts the investigation of asymptotic
properties of the metal-insulator transition, as it would require much larger
cluster sizes that are beyond computational capabilities.

Cluster Perturbation Theory, Hubbard model, metal-insulator transition

##  I Introduction

The Hubbard model probably belongs to the most studied systems in solid state
theory. Although its Hamiltonian possesses a simple form, it captures
important aspects of various many-body phenomena like Mott-insulating states,
antiferromagnetism and superconductivity.reich_heavy-fermion_1988 ;
tremblay_pseudogap_2006 ; huscroft_pseudogaps_2001 ;
gull_superconductivity_2013 ; Arovas_hubbard_model . The Hamiltonian has three
terms: the first term describes the hopping of the electrons on the lattice,
the second term a repulsive Coulomb interaction of spin up and spin down
electrons on the same site and the third term is the chemical potential, which
we shifted such that half filling corresponds to $\mu=0$ for bi-partite
lattices:

$\displaystyle{\cal H}$
$\displaystyle=-\sum_{\sigma}\sum_{x,y}t_{x,y}\hat{c}^{\dagger}_{x,\sigma}\hat{c}_{y,\sigma}+U\sum_{x}\hat{n}_{x,\uparrow}\hat{n}_{x,\downarrow}$
$\displaystyle-(\mu+U/2)\sum_{x}(\hat{n}_{x,\uparrow}+\hat{n}_{x,\downarrow}),$
(1)

where $x$ and $y$ are labelling the lattice sites and
$\sigma=\uparrow,\downarrow$ denotes the spin index. With
$\hat{c}^{\dagger}_{x,\sigma},\hat{c}_{y,\sigma}$ we denote the fermionic
creation and annihilation operators and
$\hat{n}_{x,\sigma}=\hat{c}^{\dagger}_{x,\sigma}\hat{c}_{x,\sigma}$ is the
occupation number operator.

One interesting aspect of the Hubbard model is its Mott-insulating state at
high interaction strength as well as the associated metal-insulator transition
it supposedly captures. In this regard a pseudogap regime at intermediate
interaction strength has been discussed tremblay_pseudogap_2006 ;
schafer_fate_2015 . Within this study we investigated this regime using
Cluster Perturbation Theory (CPT). Introduced by Senechal et al.
senechal_spectral_2000  CPT has shown remarkable results when applied to the
Hubbard model, despite being of low numerical cost. While these results caught
our initial interest for the method, we came to the conclusion that care has
to be taken when interpreting the results of CPT, especially concerning
features like spectral gaps. In the following, we will first outline the
method, apply it to the systems of interest and then analyse carefully the
accuracy of the results by comparing the one dimensional case to exact results
using Bethe ansatz.

##  II Methods

###  II.1 Cluster Green’s functions

In Cluster Perturbation Theory (CPT) the main objective is to construct an
approximation to the retarded Green’s function ${\cal
G}^{r}(\textbf{k},\omega)$ of a given lattice system in the thermodynamic
limit. This function is especially useful as it provides direct access to the
spectral function bruus_many-body_2004 ; fetter2012quantum :

$\displaystyle{\cal A}(\textbf{k},\omega)=-\frac{1}{\pi}\text{Im}{\cal
G}^{r}(\textbf{k},\omega).$ (2)

As CPT aims at approximating the Green’s function of the full system by
combining the solutions of small finite clusters cut out of the infinite
lattice, we first have to discuss how to obtain the interacting Green’s
function on such a cluster. For this we first define the retarded Green’s
function for two fermionic operators $\hat{A}$ and $\hat{B}$ as:

$\displaystyle{\cal{G}}^{r}_{\hat{A},\hat{B}}(t,t^{\prime})=-i\Theta(t-t^{\prime})\langle\\{\hat{A}(t),\hat{B}(t^{\prime})\\}\rangle$
(3)

where $\\{\ldots,\ldots\\}$ is the anticommutator, $t$ and $t^{\prime}$ are
time arguments and $\Theta(t)$ is the Heavyside Theta function. Note that we
are only interested in the $T=0$ case, which means that the expectation value
($\langle...\rangle$) only consists of the ground state $\ket{\Psi_{0}}$. For
the retarded Green’s function we have $t>t^{\prime}$ and as the Hamiltonian of
interest is time-independent, we can assume $t^{\prime}$ to be zero. As we are
going to use a Chebyshev expansion, it is convenient to rewrite the retarded
Green’s function in terms of two new functions
${\cal{G}}^{+}_{\hat{A},\hat{B}}(t)$ and ${\cal{G}}^{-}_{\hat{B},\hat{A}}(t)$
schmitteckert_calculating_2010 :

$\displaystyle{\cal{G}}^{r}_{\hat{A},\hat{B}}(t)$
$\displaystyle=-i\Theta(t)\langle\\{\hat{A}(t),\hat{B}(0)\\}\rangle,$
$\displaystyle=-i\Theta(t)\langle\hat{A}(t)\hat{B}(0)\rangle-i\Theta(t)\langle\hat{B}(0)\hat{A}(t)\rangle,$
$\displaystyle={\cal{G}}^{+}_{\hat{A},\hat{B}}(t)-{\cal{G}}^{-}_{\hat{B},\hat{A}}(t).$
(4)

Where we used the definitions:

$\displaystyle{\cal{G}}^{+}_{\hat{A},\hat{B}}(t)$
$\displaystyle=-i\Theta(t)\langle\hat{A}(t)\hat{B}(0)\rangle,$ (5)
$\displaystyle{\cal{G}}^{-}_{\hat{B},\hat{A}}(t)$
$\displaystyle=i\Theta(t)\langle\hat{B}(0)\hat{A}(t)\rangle.$ (6)

Performing a Fourier transformation we can obtain the Green’s function in the
frequency domain as:

$\displaystyle{\cal{G}}^{+}_{\hat{A},\hat{B}}(\omega)=-\bra{\Psi_{0}}\hat{A}[{\cal
H}-E_{0}-(\omega+i\eta)]^{-1}\hat{B}\ket{\Psi_{0}}.$ (7)
$\displaystyle{\cal{G}}^{-}_{\hat{B},\hat{A}}(\omega)=-\bra{\Psi_{0}}\hat{B}[{\cal
H}-E_{0}+(\omega+i\eta)]^{-1}\hat{A}\ket{\Psi_{0}},$ (8)

where $\eta>0$ is an infinitesimal parameter that ensures convergence.

###  II.2 Chebyshev Expansion

Expressions like (7) and (8) can be very efficiently handled using Chebyshev
polynomials braun_numerical_2014 . They contain the function

$\displaystyle f^{\pm}_{z}(x)=-i\int^{\pm\infty}_{0}{\mathrm{e}}^{i(\pm
z-x)t}\,dt=\frac{1}{\pm z-x},$ (9)

with $x,\text{Re}(z)\in\mathbb{R}$ and $\text{Im}(z)>0$, that can be expanded
using Chebyshev polynomials of the first kind $T_{n}(x)$:

$\displaystyle f^{\pm}_{z}(x)=\sum_{n=0}^{\infty}\alpha^{\pm}_{n}(z)T_{n}(x),$
(10)

with the expansion coefficients:

$\displaystyle\alpha^{\pm}_{n}(z)=\frac{2/(1+\delta_{n,0})}{(\pm
z)^{n+1}(1+\sqrt{z^{2}}\sqrt{z^{2}-1}/z^{2})^{n}\sqrt{1-1/z^{2}}}.$ (11)

For the polynomials the following recursion relation holds:

$\displaystyle\ket{\Phi_{0}}=b\ket{\Psi_{0}},$ (12)
$\displaystyle\ket{\Phi_{1}}=[a({\cal H}-E_{0})-b]\ket{\Phi_{0}},$ (13)
$\displaystyle\ket{\Phi_{n+1}}=2[a({\cal
H}-E_{0})-b]\ket{\Phi_{n}}-\ket{\Phi_{n-1}},$ (14)

where we choose the two parameters $a,b\in\mathbb{R}$ to fit the spectrum of
the operator $a({\cal H}-E_{0})-b$ into the interval $(-1,1)$, required by the
orthogonality relation of the Chebyshev polynomials. With this we can identify
the Green’s functions as:

$\displaystyle{\cal{G}}^{\pm}_{\hat{B},\hat{A}}(\omega)$
$\displaystyle=a\sum_{n=0}^{\infty}\alpha_{n}^{\pm}(\pm
a(\omega+i\eta)-b)\mu_{n},$ (15)

where the $\mu_{n}$ are often referred to as Chebyshev moments and are defined
as the expectation values of the polynomials:

$\displaystyle\mu_{n}=\bra{\Psi_{0}}\hat{A}T_{n}(a({\cal
H}-E_{0})-b)\hat{B}\ket{\Psi_{0}}=\bra{\Psi_{0}}\hat{A}\ket{\Phi_{n}}.$ (16)

In order to calculate these moments for the Green’s function of the finite
cluster, we require the Hamiltonian in a many particle basis. To construct the
ground state we employ a sparse matrix diagonalization as for example
introduced in Ref. senechal_introduction_2010 . The main idea is to explicitly
encode how a specific Hamiltonian acts on basis states in the occupation
number representation. For this, one needs to at least encode all the basis
states in a particular number sector. Although this can be efficiently done by
saving each basis state as the bitwise representation of an integer, the
computational space still grows exponentially, making it only usable for very
small clusters. Due to limited computational resources our calculations did
not exceed 18 site calculations. Having constructed the Hamiltonian the
groundstate can be calculated using a Lanczos algorithm.

###  II.3 Cluster Perturbation Theory

The goal of Cluster Perturbation Theory (CPT) is to approximate the Green’s
function of a particular lattice model in the thermodynamic limit by combining
the Green’s functions of small individual clusters, for example calculated as
described in the previous section. Introductions to this method are presented
in Ref. senechal_introduction_2010 ; pavarini_ldadmft_2011 . The first step is
to split the Hamiltonian into two parts:

$\displaystyle{\cal H}=\sum_{\alpha}{\cal H}_{\alpha}^{\text{cluster}}+{\cal
H}^{\text{inter}}.$ (17)

In the first part:

$\displaystyle{\cal H}^{\text{cluster}}_{\alpha}$ $\displaystyle=({\cal
H}^{c}_{0}+{\cal
H}^{c}_{U})_{\alpha}=-t\sum_{\sigma}\sum_{x,y\in\gamma^{c}_{\alpha}}\hat{c}^{\dagger}_{x,\sigma}\hat{c}_{y,\sigma}$
$\displaystyle+U\sum_{x\in\gamma^{c}_{\alpha}}\left(\hat{n}_{x,\uparrow}-\frac{1}{2}\right)\left(\hat{n}_{x,\downarrow}-\frac{1}{2}\right),$
(18)

one has the full Hubbard model on small, individual clusters
$\gamma^{c}_{\alpha}$, each labeled by the index $\alpha$ and in the second
part:

$\displaystyle{\cal
H}^{\text{inter}}=-t\sum_{\sigma}\sum_{x\in\gamma^{c}_{\alpha},y\in\gamma^{c}_{\beta}}\hat{c}^{\dagger}_{x,\sigma}\hat{c}_{y,\sigma},$
(19)

we only have the hopping elements between these individual clusters. Note that
due to this splitting, it can be very useful to describe any lattice site
$\textbf{R}_{i}$ by a combination of two new vectors:

$\displaystyle\textbf{R}_{i}=\textbf{r}_{\alpha}+\textbf{r}_{m},$ (20)

where $\textbf{r}_{\alpha}$ is the position of the individual clusters in a
new superlattice $\Gamma$ and $\textbf{r}_{m}$ describes the position of an
individual site within a cluster.
One can calculate the Green’s function for one of these clusters and use this
result for the Green’s function of all other clusters due to the lattice
symmetry. We will refer to this Green’s function as the cluster Green’s
function ${\cal G}^{c}(\textbf{r}_{m},\textbf{r}_{n},\omega)$. The main idea
within CPT consists in calculating the self-energy from the cluster Green’s
function and use it to construct an approximation for the self-energy of the
full system. We can obtain the cluster self-energy
$\Sigma^{c}(\textbf{r}_{m},\textbf{r}_{n},\omega)$ from a Dyson equation:

$\displaystyle\Sigma^{c}(\textbf{r}_{m},\textbf{r}_{n},\omega)$
$\displaystyle=({\cal G}_{0}^{c}(\textbf{r}_{m},\textbf{r}_{n},\omega))^{-1}$
$\displaystyle-({\cal G}^{c}(\textbf{r}_{m},\textbf{r}_{n},\omega))^{-1},$
(21)

where ${\cal G}_{0}^{c}(\textbf{r}_{m},\textbf{r}_{n},\omega)$ is the non-
interacting Green’s function on the cluster defined as:

$\displaystyle({\cal
G}_{0}^{c}(\textbf{r}_{m},\textbf{r}_{n},\omega))^{-1}=\omega+i\eta-{\cal
H}^{c}_{0}(r_{m},r_{n}).$ (22)

Therefore we obtain for a particular entry of the system self-energy
$\Sigma^{s}(\textbf{R}_{i},\textbf{R}_{j},\omega)$ in real space, connecting
two sites on the same cluster:

$\displaystyle\Sigma^{s}(\textbf{R}_{i},\textbf{R}_{j},\omega)$
$\displaystyle=\Sigma^{s}(\textbf{r}_{\alpha}+\textbf{r}_{m},\textbf{r}_{\alpha}+\textbf{r}_{n},\omega)$
$\displaystyle=\Sigma^{c}(\textbf{r}_{m},\textbf{r}_{n},\omega),$ (23)

and all entries of the self-energy connecting sites on different clusters are
set to zero.
Finally, we can use this approximation of the self energy, namely using the
cluster self energy for the self energy of the full system, in a Dyson
equation as before, to obtain the Green’s function of the full system:

$\displaystyle({\cal G}_{0}^{s}(\omega))^{-1}=\omega+i\eta-{\cal
H}^{\text{inter}}-\sum_{\alpha}{\cal H}^{c}_{0,\alpha}.$ (24)

Note that in this way, we treat the non-interacting part exactly. This is why
one should view the CPT approximation as a perturbation theory in $U$ rather
than a perturbation in the inter-cluster hopping. Finally we end up with the
following expression for the Green’s function of the full system:

$\displaystyle{\cal G}^{s}(\textbf{R}_{i},\textbf{R}_{j},\omega)$
$\displaystyle=(({\cal G}_{0}^{s}(\textbf{R}_{i},\textbf{R}_{j},\omega))^{-1}$
$\displaystyle-\Sigma^{s}(\textbf{R}_{i},\textbf{R}_{j},\omega))^{-1}$ (25)

###  II.4 Periodization

While the just described procedure works for finite systems, it is important
to note that there are also so called periodization schemes, which allow to
extend these results to infinite systems. Here we are going to use the so
called G-scheme, as discussed in Ref. avella_strongly_2012 . The main idea is
based on arranging clusters in an infinite superlattice and exploiting its
translational symmetry. As pointed out before, one can split any lattice
vector into one vector defined on the superlattice and one on a cluster.
Therefore, one can similarly split any wave vector k of the 1st BZ into a
combination of a wave vector in a reduced BZ $\tilde{\textbf{k}}$ associated
with the superlattice and one of the Brillouin zone of a single cluster K.
This also allows one to split the Fourier transformations into two parts, one
for the cluster and one for the superlattice. Using Bloch’s theorem for the
superlattice, one ends up with the following form for a periodized Green’s
function:

$\displaystyle{\cal
G}(\textbf{k},\omega)=\frac{1}{L}\sum_{a,b}e^{-i\textbf{k}(\textbf{r}_{a}-\textbf{r}_{b})}{\cal
G}_{a,b}(\tilde{\textbf{k}},\omega).$ (26)

##  III Results

Figure 1: First Brillouin zone (1. BZ) of the 2D square lattice with symmetry
points and the k-path for the bandstructure plots. Figure 2: Spectral function
of the Hubbard model with U=4t on a 2D square lattice, plotted along the high
symmetry axis of the 1. BZ. The broadening parameter was chosen as $\eta=0.5$.
The cluster calculations were performed on a 4x4 cluster with 120 Chebyshev
moments. Figure 3: Same spectral function as in Fig. 2, plotted only at the
$X$ point ($k=[\pi,0,0]$) and the midpoint between $\Gamma$ and $M$ point
($k=[\pi/2,\pi/2,0]$). We can still see a significant spectral weight at
$\omega=0$ indicating a pseudogap.

We can now use the just presented methods to calculate the spectral function
for our main system of interest, the 2D Hubbard model on a square lattice at
half filling (Fig. 1, 2, 3). As we can clearly see when considering the
spectral function at the $X$ point and in the midpoint between $\Gamma$ and
$M$, CPT here suggests a considerable spectral weight in the gap at U=4t.
However, we are going to show that there are two parameters that pose large
additional constraints on the resolution that can actually be achieved using
CPT. The first parameter is the convergence aiding factor $\eta$ and the
second is the finite cluster size. While the artefacts induced by the
convergence aiding factor are related to our approximation of using only a
finite number of Chebyshev moments, the constraints imposed by the finite
cluster size are an inherent limitation of CPT. These additional constraints
are typically not discussed in detail in the literature, but as we are going
to show they actually prohibit us from making accurate judgements about the
pseudogap at intermediate interaction strengths. To see this, we will
concentrate on the 1D Hubbard model since on the one hand it allows us to
compare to exact results from Bethe ansatz and on the other hand it gives a
higher resolution in k-space.

###  III.1 Convergence Aiding Factor

The convergence aiding factor $\eta$ enters the Green’s function in Eqs. (7)
and (8). We calculate these Green’s functions with the help of the Chebyshev
expansion (15) that, in practise, is evaluated only with a finite number of
Chebyshev moments. This leads to artefacts in the spectral function in the
form of Gibbs oscillations. This is illustrated for the spectral function of
the non-interaction 1D chain evaluated for a cluster with 16 sites, see Fig.
4. In Fig. 5 we show the spectral function for a specific k-value as a
function of frequency that clearly displays oscillations. Note that a higher
Chebyshev order only increases the frequency of these oscillations (see Fig.
6) but does not change their magnitude. These oscillations can be identified
as Gibbs oscillations that usually arise when approximating a sharp step
function (in this case the $\delta$-peak) by a finite Fourier expansion
series. In order to suppress these artificial oscillations we will choose a
sufficiently large value for the broadening parameter $\eta$, such that the
Chebyshev expansion is capable of resolving the peak without Gibbs
oscillations. In addition, the finite cluster is naturally characterized by a
finite level spacing. In order to mimic an infinite system and to obtain
smooth bands, a finite broadening parameter $\eta$ has to be chosen such that
it smears out the effect of the finite level spacing
schmitteckert_calculating_2010 .

Figure 4: Spectral function obtained from a cluster Green’s function of a 16
site tight binding chain with a broadening parameter of $10^{-7}$ leading to
negative values in the spectral function. Figure 5: Spectral function of a 16
site tight binding chain at $k=3\pi/4$ with $n_{ch}=60$ chebyshev moments and
a broadening parameter of $10^{-7}$. We can see oscillations around the peak
leading to negative values in the spectral function. Figure 6: Same spectral
function as in Fig. 4 but calculated using $n_{ch}=120$ Chebyshev moments.
Note that the amplitude of the oscillations does not decrease, but their
frequency increases in accordance with Gibbs phenomenon. Figure 7: Again the
same spectral function as in Fig. 4 calculated using $n_{ch}=120$ Chebyshev
moments but now with a broadening of $\eta=0.25$. The oscillations disappear,
because the peak is artificially broadened. Figure 8: Removing the artificial
broadening by means of eq. (27) results in a sharp $\delta$-like peak in the
spectrum.

While this means that altogether one has to choose $\eta$ rather large (Fig.
7), one should also realize that one can counteract the effects of this
broadening to a large extend by including the same large parameter for $\eta$
in the non interacting Green’s function when calculating the self energy.
Effectively this corresponds to subtracting $\eta$ from the inverse of the
Green’s function,

$\displaystyle{\cal G}^{c}(\omega)=\left[({\cal
G}^{c}(\omega,\eta_{C}))^{-1}-i\eta_{C}\right]^{-1}$ (27)

While this procedure works extremely well, as shown in Fig. 8, it still
depends on the exact choice for $\eta$, and it is a priori unclear which value
to choose for $\eta$. Here we want to propose two different approaches.

The first approach uses the typical single particle level spacing of the non
interacting cluster for the broadening parameter $\eta$, that can be estimated
as

$\displaystyle\eta=\frac{4t}{M_{C}}$ (28)

where $4t$ is the bandwidth with $t$ the hopping amplitude and $M_{C}$ the
cluster size. This smears the discrete cluster levels (Fig. 9), and leads to a
good approximation of the continuous cosine band structure one expects for a
1D system (Fig. 10). Applying this choice for the CPT approximation to the 1D
Hubbard model, one finds a finite spectral weight within the Hubbard gap
illustrating the numerical artefact that is induced by a large $\eta$, see
Figs. 11 and 12.

Figure 9: Spectral function of a 16 site Hubbard model at U=0 (tight binding
chain) with a broadening of $\eta=0.15$. We can see individual levels of the
cluster. Figure 10: Spectral function for the same system as in Fig. 9 but
with a broadening of $\eta=0.25$. We can see that the peaks start to overlap
and to resemble a continuous cosine bandstructure, as expected for an infinite
tight binding chain. Figure 11: CPT result for a 1D Hubbard model at U = 2t
based on a 16 site cluster with broadening $\eta=0.25$, based on the average
single particle level splitting. We obtain two continuous bandstructures with
the expected Hubbard gap. Figure 12: CPT result for the same system as in Fig.
11 shown only for $\mathbb{k}=\pi/2$ but with different broadening parameters.
We can see, that the residual spectral weight in the gap strongly depends on
the broadening parameter.

As a second approach we propose an extrapolation scheme that calculates the
CPT Green’s function for multiple values of $\eta$ and performs an
extrapolation of the results to $\eta=0$ (see Figs. 13 and 14).

Figure 13: Again CPT results for the same system as in Fig. 11 but using an
extrapolation scheme for the broadening parameter. We again obtain two
continuous bands separated by the Hubbard gap, however the bands are narrower.
Figure 14: Same plot as in Fig. 12 but now for the CPT data obtained with the
extrapolation scheme for the broadening parameter. We can see that the
spectral function in the gap almost disappears.

Although this procedure is physically sound and results in sharp peaks, there
is no guarantee that this procedure will lead to the correct thermodynamic
limit. In addition, obtaining a sharper peak does not automatically provide a
more accurate result. Only in cases where the actual width of the peak is
resolved by the Chebyshev expansion, i.e. wider than the many particle
bandwidth divided by the number of Chebyshev moments, the extrapolated peaks
could be considered reliable. Remarkably, despite the peak height and width
being $\eta$ dependent, the peak position appears to be rather stable against
a variation of $\eta$.

###  III.2 Accuracy of CPT results

It is important to note that the finite size level spacing of the cluster used
within CPT acts as a cutoff that limits any spectral resolution. Spectral
features like gaps can only be resolved reliably if they are larger than this
cutoff. As far as interaction effects are concerned, CPT only accurately
reflects the same information that a careful analysis of the cluster result
would also provide.

To extract such an information from the cluster result directly, let us
consider the spectral function of a 16 site 1D Hubbard chain without
interaction as was shown in Fig. 9 and with an interaction of U=4t (Fig. 15).
As one would expect for the non interacting case we can see 16 individual
levels forming a cosine shaped band, while for the interacting case the levels
in the middle of the spectrum move apart. Hence, the actual band gap is given
by the shift of the energy levels rather than their frequency difference
directly. Table 1 shows the frequency difference of a particular interaction
strength subtracted by the difference in the non interacting case, while the
table below shows the exact results one would get in the thermodynamic limit
using Bethe ansatz (Tab. 2). As we can see, only for $U=4t,8t$ the gap size
agrees up to the second decimal place with the exact result. However, this
means that for $U=1t,2t$ and using CPT with the currently computationally
accessible cluster sizes we can not judge if there is an actual gap in the
system, as the deviation is on the same order of magnitude as the actual gap.
Additionally we can see that the gap size does not change significantly with
the cluster size. Hence we have to assume that the considered cluster sizes
are simply to small to resolve the gap accurately.

Figure 15: Spectral function of a 16 site Hubbard chain at U=4t and a broadening parameter of $\eta=0.15$. Comparing to the non interacting result in Fig. 9 we can see shifts in the positions of the single particle levels. Figure 16: The average of the spectral function at the k points $k_{1/2}=\frac{M}{2}\pm 1$ for the M=16 site calculations as shown in Fig. 9.  Figure 17: The same plot as in Fig. 16 but with an interaction strength of U=4t. We can see the shift in the peak positions which corresponds to the Hubbard gap. | M | U=1t | U=2t | U=4t | U=8t
---|---|---|---|---|---
6 | $(\Delta(U)-\Delta(0))/2$ | 0.03 | 0.16 | 0.66 | 2.24
8 | $(\Delta(U)-\Delta(0))/2$ | 0.03 | 0.15 | 0.64 | 2.24
10 | $(\Delta(U)-\Delta(0))/2$ | 0.02 | 0.14 | 0.63 | 2.24
12 | $(\Delta(U)-\Delta(0))/2$ | 0.02 | 0.13 | 0.62 | 2.25
14 | $(\Delta(U)-\Delta(0))/2$ | 0.02 | 0.13 | 0.62 | 2.26
16 | $(\Delta(U)-\Delta(0))/2$ | 0.03 | 0.13 | 0.63 | 2.30
Table 1: Half of the band gap $\Delta$, calculated using the peak positions
from the interacting case and correcting them using the peak positions of the
non interacting one.

Now in order to arrive at the same result using the full CPT calculation, we
can simply plot the spectral function in the middle of the spectrum. If we do
this for multiple cluster sizes (Fig. 18), we can see, that the peaks slightly
move to the center with increasing cluster size. Hence we can extrapolate the
peak positions for the gap of an infinite cluster as shown in Fig. 19. This
results in gaps comparable to results obtained directly from the cluster (Tab.
1). As in the case where we just considered the cluster we see that the error
is too large to judge the gap accurately at $U=1t,2t$.

Figure 18: Shown is the influence of the cluster size on the band gap for $U/t=4$. We show the spectral function at $k=\pi/2$. The peaks defining the gap get broader and move closer together the larger the cluster. Note that the minimum spectral weight within the gap stays almost constant with system size and is almost zero due to our choice of $\eta=0.15$. The peaks which are marked are used for a finite size extrapolation in Fig. 19. Figure 19: The $\omega$ values of the peaks marked in Fig. 18 are plotted against the inverse of the corresponding cluster size. Doing so allows for a finite size extrapolation to zero, corresponding to an infinite cluster size. The arrangement of the data suggests a linear extrapolation, which we performed and the result at zero is shown in the title. It corresponds to half the bandgap ($\Delta_{CPT}/2$). The second value was obtained, by doing the same extrapolation for the deviation from the exact Bethe result ($\Delta_{CPT}/2-\Delta_{Bethe}/2$). | U=1t | U=2t | U=4t | U=8t
---|---|---|---|---
$\Delta_{\text{Bethe}}/2$ | 0.003 | 0.086 | 0.643 | 2.340
$\Delta_{CPT}/2$ | 0.028 | 0.143 | 0.683 | 2.344
$(\Delta_{CPT}/2)_{Error}$ | 0.025 | 0.057 | 0.040 | 0.004
Table 2: In row one, we show half the band gap $\Delta$ calculated using Bethe
ansatz. In row two, we give half the gap size obtained via CPT and in row
three the extrapolation of the deviation. We can see that the results deviate
on the order of magnitude of $10^{-1}$, as they already did, when we were just
considering the cluster result. This was expected, since the accuracy of CPT
is determined by the finite cluster.

After this discussion, we now return to the 2D case which sparked this
investigation in the first place. We come to the conclusion that CPT is not
able to resolve reliably the spectral weight within the gap and, in
particular, to predict the presence or absence of a pseudogap despite the
resolution that the plot in Fig. 2 suggests. This is, firstly, due to the
influence of the finite broadening parameter used to dampen the Gibbs
oscillation induced by the Chebyshev approximation, and, secondly, due to the
finite level spacing associated with the chosen 2D cluster.

###  III.3 Use Cases for CPT

One might now pose the question as to why one should employ CPT in the first
place, since all of the reliable information is already contained in the
cluster results. To this end, one should realize that just judging from the
cluster results it can be very hard to identify how the actual band structure
of a system would look like. CPT though, acting as a sort of interpolation
scheme between the single particle levels of the cluster and the infinitely
large lattice, can give a very good first guess of the band structure, which
comes at a negligible additional computational cost. Especially for materials,
where short range correlations are predominant, this guess will be highly
accurate. Hence one might view it as a useful tool when scanning through a
variety of materials. One can simply perform fast CPT calculations for each
material, identify interesting band structures and then use more advanced
methods to investigate these materials further. As an example, we provide the
spectral functions for a more exotic lattice, the Lieb lattice
weeks_topological_2010  in Fig. 20 and 21.

Figure 20: Spectral function of the Lieb lattice at $U=0$ obtained via CPT
using a 2x2 cluster of three atomic unitcells. We can see the flat band at
$\omega=0$. Figure 21: The CPT result for the spectral function of the Lieb
lattice, but at $U=4t$. We can observe a splitting of the flat band.

##  IV Conclusion

In this study we have presented fundamental limitations of the CPT method.
Most importantly we showed that the cluster Green’s function already contains,
as far as interaction effects are concerned, all the information that is
included in the CPT Green’s function and CPT only stays consistent with this
information. Analysing the 1D case we have shown that the current
computational limitations prohibit us to make accurate judgements of features
like the Hubbard gap at intermediate interaction strengths, as this would
require calculations on larger clusters, as the resolution is limited by the
finite size induced level splitting. Additionally and specifically, for the
approach of calculating the cluster Green’s function via a Chebyshev
expansion, we have shown that Gibbs oscillations require the choice of a large
broadening parameter that in return prohibits one from making accurate
judgements about the broadness of the peaks in the single particle spectrum.
Hence we conclude that use cases for CPT are limited to cases where one is
interested in obtaining numerically cheap initial guesses for the spectral
function of materials with short range correlations, while more advanced
methods need to be employed, to gain higher resolution.

Data availability The data presented in this publication is avaliable on
Zenodo under the DOI: 10.5281/zenodo.8063247.

Funding This project was supported by BMBF via the MANIQU grant no. 13N15576.

## References

  * [1] Ariel Reich and L. M. Falicov.  Heavy-fermion system: An exact many-body solution to a periodic-cluster hubbard model.  Phys. Rev. B, 37(10):5560–5570, 1988.
  * [2] A.-M. S. Tremblay, B. Kyung, and D. Sénéchal.  Pseudogap and high-temperature superconductivity from weak to strong coupling. towards quantitative theory.  Low Temperature Physics, 32(4):424–451, 2006.
  * [3] C. Huscroft, M. Jarrell, Th. Maier, S. Moukouri, and A. N. Tahvildarzadeh.  Pseudogaps in the 2d hubbard model.  Phys. Rev. Lett., 86(1):139–142, 2001.
  * [4] Emanuel Gull, Olivier Parcollet, and Andrew J. Millis.  Superconductivity and the pseudogap in the two-dimensional hubbard model.  Phys. Rev. Lett., 110(21):216405, 2013.
  * [5] Daniel P. Arovas, Erez Berg, Steven A. Kivelson, and Srinivas Raghu.  The hubbard model.  Annual Review of Condensed Matter Physics, 13(1):239–274, 2022\.
  * [6] T. Schäfer, F. Geles, D. Rost, G. Rohringer, E. Arrigoni, K. Held, N. Blümer, M. Aichhorn, and A. Toschi.  Fate of the false mott-hubbard transition in two dimensions.  Phys. Rev. B, 91(12):125109, 2015.
  * [7] D. Sénéchal, D. Perez, and M. Pioro-Ladrière.  Spectral weight of the hubbard model through cluster perturbation theory.  Phys. Rev. Lett., 84(3):522–525, 2000.
  * [8] Henrik Bruus and Karsten Flensberg.  Many-Body Quantum Theory in Condensed Matter Physics: An Introduction.  OUP Oxford, 2004.
  * [9] Alexander L Fetter and John Dirk Walecka.  Quantum theory of many-particle systems.  Courier Corporation, 2012.
  * [10] Peter Schmitteckert.  Calculating green functions from finite systems.  J. Phys.: Conf. Ser., 220:012022, 2010.
  * [11] A. Braun and P. Schmitteckert.  Numerical evaluation of green’s functions based on the chebyshev expansion.  Phys. Rev. B, 90(16):165112, 2014.
  * [12] David Sénéchal.  An introduction to quantum cluster methods.  arXiv:0806.2690 [cond-mat], 2010.
  * [13] Eva Pavarini, Erik Koch, Dieter Vollhardt, Alexander I. Lichtenstein, Institute for Advanced Simulation, German Research School for Simulation Sciences, and Deutsche Forschungsgemeinschaft, editors.  The LDA+DMFT approach to strongly correlated materials: lecture notes of the Autumn School 2011 Hands-on LDA+DMFT.  Number Band 1. 2011.
  * [14] Adolfo Avella and Ferdinando Mancini, editors.  Strongly Correlated Systems: Theoretical Methods, volume 171 of Springer Series in Solid-State Sciences.  Springer Berlin Heidelberg, 2012.
  * [15] C. Weeks and M. Franz.  Topological insulators on the Lieb and perovskite lattices.  Physical Review B, 82(8):085310, August 2010.  Publisher: American Physical Society.