Accuracy of the recovered flux of extended sources obscured by bad pixels in
the central EPIC FOV

XMM-SOC-CAL-TN-0227

J. Nevalainen, I. Valtchahov, R. D. Saxton and S. Molendi

Feb 2 2021

##  1 Introduction

A fraction of the XMM-Newton EPIC FOV is obscured by the dysfunctional (i.e.
bad) pixels, dominated by the CCD gaps, which complicates the analysis of
extended X-ray sources. Before SAS version 17.0.0 the only standard option for
attempting the recovery of the obscured flux was to scale the measured flux by
the fraction of the obscured detector area. This is accurate only in the very
limiting case of an extended source with a uniform spatial flux distribution.

Since version 17.0.0 SAS has included an option for recovering the lost
fraction of the flux obtained by a primary instrument (e.g. pn) by utilising
the information of the spatial distribution of the flux (i.e. an image) of the
same source observed with another instrument (e.g. MOS2). Assuming that the
supplementary image is minimally obscured at the spatial locations where the
primary instrument is obscured, the correction can be accurate.

The procedure is implemented to arfgen whereby the effective area associated
with the input spectrum is reduced by the lost flux fraction in the intended
full extraction region, as informed by the supplementary image. While it is
the flux, rather than the effective area, that is reduced in reality, the
above correction works in practice. Namely, when using the effective area,
reduced by the above procedure, in the consequent spectral fitting of the data
reduced by the bad pixel obscuration by the same factor, the net effect is the
increase of the spectral emission model normalisation compared to that
obtained with the original effective area. The increased model normalisation
corresponds to the full unobscured extraction region, as desired.

Since the fraction of the obscured area may exceed 10% in the EPIC FOV (see
below) the knowledge of the accuracy of the recovery procedure is important if
one is interested in a high level of accuracy of the recovered flux of an
extended source in the full intended extraction region. Our aim in this note
is to evaluate the accuracy achieved by the recommended default procedure
using a sample of clusters of galaxies. We do not examine here the accuracy of
the recovery of the spectral shape since its change is expected to be small in
general compared to the flux effect.

##  2 Data

###  2.1 The sample

We used a sample of 27 XMM-Newton observations of nearby (z $<$ 0.1) clusters
of galaxies for our work (see Table 1). Galaxy clusters are useful because
their stable, bright and diffuse X-ray emission fills completely the central
regions of the EPIC FOV, providing sufficient numbers of counts within the
pixels of our adopted size of 5 arcsec (see Section 3.1). Also, the surface
brightness of clusters decreases quite rapidly with the distance from the
cluster center (see Section 3.2.2) which is essential for the testing of the
recovery of the obscured non-uniform flux. We limit the study within 6 arcmin
distance from the cluster center in order to maintain the background below 10%
of the cluster emission in the adopted 0.5-7.0 keV energy band.

###  2.2 Processing

When processing the data and extracting products for generating the effective
area files (i.e. event files, spectra and images), we intended to mimic the
choices of a general user for wider applicability of our results, i.e. we use
default procedures and parameter values. We processed the pn and MOS event
files with epchain and emchain procedures, respectively, in SAS18.0.0 with
default parameters. Relevant for our analysis, the epchain default setting
runbadpixfind=Y will cause the badpixfind routine to be executed. In the
default mode 1 (enabled by searchbadpix=Y) hot, dead and flickering pixels in
a given observation are identified. In case of emchain, the default setting
(badpixfindalgo=EM) activates a more sensitive embadpixfind routine, instead
of badpixfind routine used for pn.

The default settings runbadpix=Y for epchain and withbadpix=yes for emchain
activate the badpix routine with getnewbadpix=Y whereby the bad pixels
identified above will be propagated to the BADPIX extension of the event
files. The default settings getuplnkbadpix=Y and getotherbadpix=Y propagate
the uplinked and non-uplinked bad pixels to the BADPIX extension.

We then filtered the pn and MOS event files with the recommended expressions

    FLAG==0 && PATTERN<=4

and

    #XMMEA_EM && PATTERN<=12

We further filtered the event files by removing the data obtained during flare
periods. These event files (e.g. pn.fits) will be used as arguments for the
badpixlocation parameter when running arfgen in Section 3.1.

We extracted spectra with default parameters within central r = 6 arcmin
region (defined in sky coordinates) centered at the X-ray peak (see Table 1).
We did not exclude point sources except in the case of A2142, where an
exceptionally bright source is partially obscured by a pn CCD gap. These
spectra (e.g. pn.pha) will be used as arguments for the spectrumset parameter
when running arfgen in Section 3.1.

We extracted raw count MOS2 images within central r = 6 arcmin region in the
sky coordinate system with a bin size of 5 arcsec in the 0.5-7.0 keV band.
These images (e.g. MOS2.im, see Fig. 1) will be used as arguments for the
detmaparray parameter when running arfgen in Section 3.1.

In order to obtain estimates of the geometric fraction Rarea of the unobscured
area to that of the full intended circle with r = 6 arcmin, we utilised the
backscale \- routine with withbadpixcorr=yes, badpixelresolution=1 (1 arcsec)
and badpixlocation set to the event file produced above (e.g. pn.fits). The
computations yielded that on average the pure geometric area reduction factor
due to bad pixel/CCD gap obscuration (1 -$<$R${}_{area}>$) is 14% for pn, 5%
for MOS1 and 4% for MOS2 for our sample (see Table 1).

Figure 1: Raw count images of A2029 cluster (OBSID = 0551780301) using pn
(left panel), MOS1 (middle panel) and MOS2 (right panel) in sky coordinates in
the central r = 6 arcmin region. The color bars indicate the level of counts
per pixel.  Table 1: Results

cluster | Obs. ID | RA | DEC |  Rareaa |  Rimageb |  Rmodelc
---|---|---|---|---|---|---
|  |  |  | ———————— |  |
|  |  |  | pn | M1 | M2 |  |
A85 | 0723802101 | 10.460333 | -9.303008 | 0.866 | 0.932 | 0.963 | 0.8819 | 0.8851
A85 | 0723802201 | 10.460333 | -9.303008 | 0.866 | 0.933 | 0.965 | 0.8807 | 0.8846
A401 | 0112260301 | 44.740833 | 13.582642 | 0.862 | 0.960 | 0.966 | 0.8564 | 0.8496
A478 | 0109880101 | 63.354875 | 10.465269 | 0.860 | 0.965 | 0.972 | 0.8828 | 0.8898
A644 | 0744412201 | 124.356792 | -7.512367 | 0.867 | 0.931 | 0.960 | 0.8962 | 0.8837
A754 | 0136740101 | 137.333375 | -9.679144 | 0.849 | 0.969 | 0.972 | 0.8607 | 0.8595
A1650 | 0093200101 | 194.673000 | -1.761561 | 0.847 | 0.958 | 0.958 | 0.8327 | 0.8287
A1651 | 0203020101 | 194.843583 | -4.196039 | 0.867 | 0.960 | 0.965 | 0.8758 | 0.8811
A1795 | 0097820101 | 207.220833 | 26.590278 | 0.887 | 0.962 | 0.965 | 0.9087 | 0.9035
A2029 | 0551780301 | 227.734167 | 5.744639 | 0.871 | 0.899 | 0.963 | 0.9312 | 0.9278
A2142 | 0674560201 | 239.583208 | 27.233158 | 0.865 | 0.945 | 0.951 | 0.8659 | 0.8742
A2244 | 0740900101 | 255.676958 | 34.060236 | 0.866 | 0.929 | 0.961 | 0.8906 | 0.8917
A2319 | 0302150101 | 290.297333 | 43.949839 | 0.863 | 0.960 | 0.960 | 0.8752 | 0.8779
A2319 | 0302150201 | 290.297333 | 43.949839 | 0.873 | 0.859 | 0.965 | 0.8901 | 0.8895
A2319 | 0600040101 | 290.297333 | 43.949839 | 0.866 | 0.954 | 0.962 | 0.8698 | 0.8643
A3112 | 0603050101 | 49.490417 | -44.238417 | 0.866 | 0.947 | 0.960 | 0.9070 | 0.9090
A3391 | 0505210401 | 96.585067 | -53.693201 | 0.866 | 0.959 | 0.967 | 0.8511 | 0.8622
A3558 | 0107260101 | 201.987125 | -31.495583 | 0.868 | 0.966 | 0.966 | 0.8750 | 0.8788
A3571 | 0086950201 | 206.867667 | -32.864533 | 0.867 | 0.965 | 0.968 | 0.8745 | 0.8713
A3827 | 0149670101 | 330.469708 | -59.946483 | 0.865 | 0.962 | 0.969 | 0.8723 | 0.8695
Coma | 0300530101 | 194.954625 | 27.960983 | 0.858 | 0.906 | 0.960 | 0.8686 | 0.8590
Coma | 0153750101 | 194.954625 | 27.960983 | 0.864 | 0.968 | 0.973 | 0.8661 | 0.8665
CygA | 0302800101 | 299.868042 | 40.733919 | 0.867 | 0.954 | 0.962 | 0.9073 | 0.9022
PKS0745 | 0105870101 | 116.880417 | -19.294394 | 0.865 | 0.961 | 0.965 | 0.9055 | 0.9146
Ophiuchus | 0505150101 | 258.115500 | -23.368714 | 0.865 | 0.955 | 0.961 | 0.8671 | 0.8630
Triangulum | 0093620101 | 249.583167 | -64.356086 | 0.854 | 0.964 | 0.970 | 0.8525 | 0.8469
ZwCl1215 | 0300211401 | 184.421333 | 3.655767 | 0.866 | 0.944 | 0.962 | 0.8637 | 0.8696

a The ratio of the areas of the good unobscured region as reported by
backscale - routine in SAS 18.0.0 to the full geometric extraction region (a
circle with r = 6 arcmin).
b The ratio of the effective area related to the pn spectrum reduced by arfgen
when using the MOS2 image for correcting for the bad pixels to that obtained
when applying no correction for the bad pixels (Eq. 1).
c The expected value of the effective area re-scaling factor, obtained via
modelling the cluster image (Eq. 3, i.e. arfgen is not involved).

##  3 Testing

Since the obscured fraction of the central r = 6 arcmin area is larger for pn
than for MOS (see Table 1 and Fig. 1), the inaccuracy of the total recovered
flux is probably largest for pn. We performed our test only for the pn flux
recovery to get a single, conservative estimate for the accuracy of the
procedure for all EPIC instruments.

###  3.1 Re-scaling factors (Rimage)

The relative accuracy of the recovered total flux is given by the accuracy of
the re-scaling factors (Rimage). Our goal in this work is to evaluate the
accuracy of Rimage estimated by arfgen for recovering the obscured pn flux
using MOS2 images. For this we first produced the effective area files using
syntax:

    arfgen spectrumset=pn.pha arfset=pn.arf detmaptype=dataset
    detmaparray=MOS2.im withbadpixcorr=yes badpixlocation=pn.fits
    badpixmaptype=dataset extendedsource=yes

where pn.pha, pn.fits and MOS2.im are the pn spectrum, pn event file and MOS2
image, respectively, produced as described in Section 2.2. We applied the
default bad pixel resolution of 2 arcsec. We used the resulting effective area
files (pn.arf) for each observation to obtain the re-scaled effective areas
Aeff,image as a function of the photon energy.

In order to evaluate the effective area re-scaling factors estimated and
implemented by arfgen in the above run, we also generated the effective areas
Aeff,ref by applying no bad pixel correction (badpixcorr=no) for each
observation. We then computed the re-scaling factors

$R_{image}=A_{eff,image}/A_{eff,ref}$ (1)

as a function of photon energy. In our adopted 0.5–7.0 keV band, the ratios
are almost constant (the standard deviation $\sim$0.1%). We adopted the mean
Rimage value in this band as the re-scaling factor for a given observation
(see Table 1).

On average the pn flux is reduced by $\sim$12% due to the bad pixel/CCD gap
obscuration in our sample, as estimated by arfgen, using MOS2 image for the
correction. This is a bit smaller than the purely geometric area obscuration
fraction (14%, see Section 2.2). This is due to the radially decreasing
surface brightness of galaxy clusters which renders the loss of the detector
area at larger radii less important for the total flux.

###  3.2 The “true” values (Rmodel)

We need to compare the effective area ratios Rimage obtained above to the
”true” values (Rmodel) so that we can test the accuracy of the procedure. We
obtained these by generating 1) a model pn image of a given cluster so that we
know a priori the unobscured count rate at each pixel in the full intended
extraction region and 2) the spatial mask of pn bad pixels so that we know
which pixels of the model image we should reject, in order to model the bad
pixel obscuration effect.

####  3.2.1 The masks

When modelling the bad pixel obscuration we must make a choice of the pixel
size of our masks and model images. In order to accurately 1) map the bad
pixels and 2) match the spectrum extraction region (Section 2.2) we adopted
the smallest useful pixel size of (0.05 arcsec)2, which is the size of a
single pixel in the internal XMM-Newton detector focal plane coordinate
system. This oversamples the physical EPIC-pn pixel by a factor of
82$\times$82.

We utilised the eexpmap routine of SAS for creating the spatial masks of the
bad pixels. While we are not interested in the exposure maps as such (which
this tool is designed to produce), this tool anyhow provides a practical way
of carrying out our task in hand. The idea is that bad pixels are traced by
zero values in the exposure maps. In order to avoid possible complications
arising when projecting the data from the detector plane (where the data are
observed) to the more convenient sky coordinate system preferred by general
users, we produced our bad pixel masks in the detector coordinate system.

We thus produced pn images in detector coordinates with a pixel size of 0.05
arcsec and used those as arguments for the imageset parameter when running
eexpmap. We used attrebin=0.021 (arcsec) to make sure that if the attitude
shifts by more than half the pixel size of 0.05 arcsec the translation between
detector and sky coordinates is recalculated. We turned the resulting exposure
maps into bad pixel masks by setting all pixels with zero exposure time as
zero, while all other pixel values were set to unity (see Fig. 2).

####  3.2.2 The model images

We base the synthetic pn model images on the surface brightness profile models
obtained by fitting the MOS2 imaging data of a given cluster. We divided the
raw count images by the exposure map to obtain vignetting-corrected count rate
maps. Using those we accumulated the total count rate in co-centric annuli
centered at the same point (the X-ray peak) as the spectrum extraction region
(see Section 2.2) except that we ignored the bad pixels and brightest point
sources. We divided the count rates with the area covered by the used pixels
in a given annulus to obtain surface brightness profile. The background is
defined as the constant surface brightness level at larger radii. We fitted
the profiles with a single-$\beta$ or double-$\beta$ model

$I(r)=I_{0,1}\times\left[1+{\left(\frac{r}{r_{\mathrm{core,1}}}\right)}^{2}\right]^{(-3\beta+\frac{1}{2})}+I_{0,2}\times\left[1+{\left(\frac{r}{r_{\mathrm{core,2}}}\right)}^{2}\right]^{(-3\beta+\frac{1}{2})}+bkg,$
(2)

where r is the projected radius and bkg is the constant background (see Table
2 and Fig. 2)

Figure 2: Left panel: The MOS2 surface brightness profile data (black
diamonds) of A2029 and the best-fit model (blue line) consisting of two
$\beta$-model components (red and green lines) and a constant background
(black dotted line). Middle panel: The simulated model image of A2029 using
the best-fit surface brightness model within the central r = 6 arcmin region.
Right panel: the mask used for filtering the simulated image with bad pixels
(black).  Table 2: Best-fit surface brightness profile model parameters

cluster |  I0,1 |  rcore,1 |  I0,2 |  rcore,2 | $\beta$
---|---|---|---|---|---
|  c s-1 arcmin-2 | arcmin |  c s-1 arcmin-2 | arcmin |
A85 | 0.32 | 1.6 | 1.7 | 0.31 | 0.53
A401 | 0.09 | 4.4 | 0.16 | 2.0 | 0.89
A478 | 0.88 | 1.2 | 2.3 | 0.30 | 0.65
A644 | 0.11 | 3.3 | 0.35 | 1.4 | 0.81
A754 | 0.02 | 11.2 | 0.13 | 2.4 | 0.64
A1650 | 0.13 | 2.9 | 0.52 | 0.82 | 0.91
A1651 | 0.14 | 2.9 | 0.31 | 1.1 | 0.90
A1795 | 0.27 | 1.7 | 2.0 | 0.54 | 0.59
A2029 | 1.3 | 0.95 | 2.3 | 0.27 | 0.58
A2142 | 0.43 | 1.2 | 0.55 | 0.60 | 0.53
A2244 | 0.40 | 1.2 | 0.28 | 0.33 | 0.66
A2319 | 0.15 | 2.7 | 0.23 | 1.2 | 0.49
A3112 | 0.48 | 1.1 | 2.5 | 0.24 | 0.65
A3391 | 0.04 | 4.4 | 0.04 | 1.1 | 0.74
A3558 | 0.07 | 7.3 | 0.19 | 1.6 | 1.1
A3571 | 0.26 | 3.1 | 0.28 | 0.77 | 0.60
A3827 | 0.07 | 4.0 | 0.20 | 1.6 | 1.0
Coma | 0.14 | 7.5 | — | — | 0.51
CygA | 3.5 | 0.30 | — | — | 0.51
Ophiuchus | 0.47 | 4.1 | 1.0 | 0.46 | 0.59
PKS0745 | 0.49 | 1.2 | 3.2 | 0.39 | 0.68
Triangulum | 0.06 | 8.4 | 0.23 | 2.7 | 0.77
ZwCl1215 | 0.07 | 3.8 | 0.05 | 1.5 | 0.82

The best-fit parameters for the $\beta$ model (Eq. 2) obtained by fitting the
MOS2 data.

We then simulated an unobscured image of a given cluster using the surface
brightness profile model obtained above (see Fig. 2). As in the case of the
bad pixel masks (see Section 3.2.1) we produced the cluster images in the
detector coordinate system using a pixel size of 0.05 arcsec. The synthetic
image has identical attitude and offset between the cluster and FOV centers to
those of the observational pn image and the mask image. The synthetic image
yields the model count rate $CR(model,full)$ from the unobscured full central
r=6 arcmin region of a given cluster.

The real galaxy clusters have a varying degree of deviations from the
azimuthal symmetry we assumed when modelling the clusters. We essentially
interpolate the surface brightness model to estimate the lost flux at the
locations of the bad pixels. Given the random nature of the deviations from
the symmetry we do not expect a significant bias in our modelling. This may
cause some scatter between the modelled and arfgen \- estimated fluxes which
will manifest itself in our results.

###  3.3 Results

We then multiplied the model images with the masks in order to model the
effect of the bad pixel/CCD gap obscuration. Using the mask-filtered images we
then obtained the model count rate $CR(model,unobsc)$ in the surviving
unobscured region. We then calculated the “true” flux reduction factor, i.e.
the “true” effective area re-scaling factor

$R_{model}=CR(model,unobsc)/CR(model,full)$ (3)

for each observation. We then compared these with the effective area re-
scaling factors Rimage (Eq. 1) produced with arfgen using the observational
data (see Section 3.1). The comparison yielded the following results for the
central r = 6 arcmin regions of the EPIC FOV using SAS 18.0.0. (see Fig. 3 and
Table 1).

  * •

On average, the flux reduction fraction as given by 1) the simulations (in
detector coordinate system with a pixel size of 0.05 arcsec) and by 2) arfgen
on the observational data (in sky coordinates with bad pixel resolution of 2
arcsec) differ by $\approx$0.03%

  * •

In some individual cases the difference exceeds 1% (1.4% at the maximum). The
standard deviation of the full sample is $\approx$0.7%.

  * •

The above results do not change significantly when improving the bad pixel
resolution from the default 2 arcsec to 0.2 arcsec when running arfgen.
Therefore we recommend to use the default value of 2 arcsec when analysing
extended sources.

Figure 3: The ratios of the test values (Rimage) to the model values (Rmodel)
for each observation (blue diamonds), together with the sample average
(denoted with black dotted line). Unity corresponds to a perfect match. The
red dashed lines indicate the standard deviation.

##  4 Conclusions

We examined the accuracy of the recovery method of the EPIC flux obscured by
the bad pixels in the central r = 6 arcmin region as implemented in SAS 18.0.0
i.e. using a supplementary raw count image for flux information. We used a
sample of galaxy clusters for testing. We found that the accuracy of the
recovered total flux is better than 0.1% on average while in individual cases
the recovered flux may be uncertain by $\sim$1%.

Acknowledgements
The work is supported by the Horizon 2020 Framework Programme of the European
Union Grant Agreement No. 871158. We acknowledge the support by the Estonian
Research Council grants IUT40-2, PRG1006, and by the European Regional
Development Fund (TK133). J.N. acknowledges ESAC Science Faculty support for a
short visit when this work was initiated.


# Direct Geometric Probe of Singularities in Band Structure

Charles D. Brown1,2,4∗, Shao-Wen Chang1,2, Malte N. Schwarz1,2,
Tsz-Him Leung1,2, Vladyslav Kozii1,3, Alexander Avdoshkin1,
Joel E. Moore1,3, Dan Stamper-Kurn1,2,3∗

1Department of Physics, University of California, Berkeley
366 Physics North, Berkeley, CA 94720, USA

2Challenge Institute for Quantum Computation
University of California, Berkeley CA 94720, USA

3Materials Sciences Division
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA

4Department of Physics, Yale University
217 Prospect Street, New Haven, CT 06520, USA

∗To whom correspondence should be addressed;
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
> A quantum system’s energy landscape may have points where multiple energy
> surfaces are degenerate, and that exhibit singular geometry of the
> wavefunction manifold, with important consequences for the system’s
> properties. Ultracold atoms in optical lattices have been used to indirectly
> characterize such points in the band structure. Here, we measure the non-
> Abelian transformation produced by transport directly through the
> singularities. We accelerate atoms along a quasi-momentum trajectory that
> enters, turns, and then exits the singularities at linear and quadratic band
> touching points of a honeycomb lattice. Measurements after transport
> identify the topological winding numbers of these singularities to be 1 and
> 2, respectively. Our work introduces a distinct method for probing
> singularities that enables the study of non-Dirac singularities in
> ultracold-atom quantum simulators.

## Main Text:

Energy surfaces are used to describe the structure and dynamics of quantum
systems whose Hamiltonians contain one or more continuous parameters. Notable
examples include band structure, which describes the motion of single
particles within a crystal as a function of their quasi-momentum, and the
potential energy surfaces that describe molecules as a function of their
nuclear coordinates. Each point on an energy surface corresponds to an
eigenenergy and an eigenstate of the physical system. Although the energies
themselves are highly important for explaining material (?, ?) and chemical
(?, ?, ?, ?) properties, so too are the local geometry and global topology of
the eigenstate manifolds.

The geometry of an eigenstate manifold can be revealed through transport of a
quantum state along a smooth path of parameters that define the system’s
Hamiltonian. This transport is generally nonholonomic, meaning that the state
generated by transport from an initial to a final point depends on the path
along which the system was transported. Such transport has been explored
mainly in the two limiting cases in which the energy spectrum of a system is
either largely gapped (?, ?) or entirely gapless (?) along a closed loop in
parameter space. In the former limit, the state-space geometry generates a
Berry phase; in the latter, the nonholonomy generalizes to a Wilson loop
operator describing a path-dependent rotation within the degenerate subspace.
In terms of $|u^{m}_{{\bf q}}\rangle$, the cell-periodic of the Bloch
wavefunction that describes a single particle, both types of dynamics derive
from the Berry connection matrix, ${\bf A}^{nm}_{{\bf q}}\equiv i\langle
u^{n}_{{\bf q}}|\partial_{{\bf q}}|u^{m}_{{\bf q}}\rangle$, which expresses
the local geometry of state space. Here, focusing on the case of band
structure, $n$ and $m$ are band indices, and $\mathbf{q}$ is the quasi-
momentum. In the gapped limit, the Berry phase is determined solely by one
(Abelian) diagonal element of this matrix; in the gapless limit, off-diagonal
elements enter, leading to non-Abelian state rotations (?).

In this work, we explore the nonholonomy of transport through a state space
containing singular points of degeneracy. One example of such points, which we
probe experimentally, is the Dirac points of degeneracy between the $n=1$ and
$n=2$ bands of the two-dimensional honeycomb lattice, lying at the
$\mathbf{K}$ and $\mathbf{K}^{\prime}$ points of the Brillouin zone (Fig. 1).
Away from these points, the energy gap between the touching bands grows
linearly with quasi-momentum. The singular state geometry around each linear
band touching point (LBTP) has profound implications for the material
properties of graphene, e.g. related to Klein tunneling of electrons through
potential barriers (?) and the appearance of a half-integer quantum Hall
effect (?). The Dirac point of the honeycomb lattice has been explored also in
ultracold-atom experiments (?), including interferometric measurements of the
Berry phase produced along trajectories that circle the Dirac point (?) and
direct mapping of the Bloch-state structure across the Brillouin zone (?, ?,
?).

Crystalline materials may also host a singular quadratic band-touching point
(QBTP), about which the energy gap between two bands grows quadratically with
quasi-momentum. As before, the QBTP can profoundly affect material properties.
For example, the singular QBTP is predicted to produce an anomalous Landau
level spectrum (?). Interactions can destabilize a QBTP, leading to
topologically-protected edge states, nematic phases and both quantum anomalous
Hall and spin phases (?, ?, ?, ?). The role of QBTPs is being investigated
intensely in both untwisted and twisted bilayer graphene (?, ?, ?, ?, ?).
Despite their importance, QBTPs have remained unexplored in ultracold-atom
systems.

Here, using ultracold atoms within an optical lattice, we demonstrate that
transport of a quantum state through a singular band touching point leads to a
non-Abelian, coherent state rotation between bands, with the rotation
depending on the relative orientation of path tangents entering and exiting
the singular point. Further, we show that this dependence characterizes and
distinguishes the Bloch-state geometry surrounding linear and quadratic band-
touching singularities.

Figure 1: Experimental scheme. (A) Illustration of an optical honeycomb
lattice with two sites (A and B) in the unit cell, formed by overlapping three
$\lambda=1064$ nm wavelength light beams (red arrows). Offsetting the optical
frequencies of two lattice beams by $\delta\omega_{1,2}(t)$ accelerates the
lattice and drives lattice-trapped atoms through a trajectory in quasi-
momentum. (B) The $n=\\{1,2,3,4\\}$ Brillouin zones of the honeycomb lattice
are shown in green, blue, red and purple, respectively. (C) The band structure
of the honeycomb lattice (plotted with potential depth of
$20~{}\mathrm{kHz}\times h$) exhibits an LBTP in the $s-$orbital band manifold
at $\mathbf{q}=\mathbf{K}$, and a QBTP in the $p-$orbital band manifold at
$\mathbf{q}=\mathbf{\Gamma}$.

Let us exemplify our approach by considering the $s$-orbital LBTP of a two-
dimensional honeycomb optical lattice (Fig. 1). To probe this Dirac point, we
prepare an optically trapped 87Rb Bose-Einstein condensate and then slowly
ramp up an overlain static honeycomb lattice, placing the atoms initially at
the $\mathbf{\Gamma}$-point of the $n=1$ band. Next, we apply a fictitious
force to the gas by accelerating the optical lattice potential to a velocity
$\mathbf{v}_{\mathrm{lat}}(t)$. Although the atoms remain at zero quasi-
momentum in the laboratory frame, they evolve to non-zero velocity
$\mathbf{v}=\hbar\mathbf{q}/m=-\mathbf{v}_{\mathrm{lat}}$ in the lattice
frame, with $\mathbf{q}$ being the lattice-frame quasi-momentum and $m$ being
the atomic mass.

To demonstrate the nonholonomy generated by the LBTP, we accelerate the atoms
on a trajectory that proceeds at constant acceleration from $\mathbf{\Gamma}$
to $\mathbf{K}$ (at quasi-momentum $\mathbf{q}_{\mathbf{K}}$), and thence at a
different constant acceleration to 24 equally-spaced points on a circle that
lie at distance of $0.4~{}\norm{\mathbf{q}_{K}}$ from the $\mathbf{K}$-point.
The turning angle between the rays entering and leaving the
$\mathbf{K}$-point, as defined in Fig. 2A, is varied over $\theta\in[0,2\pi]$.
We then perform “band mapping” by smoothly ramping off the lattice potential
at the fixed final quasi-momentum, i.e. with the lattice at a constant final
laboratory-frame velocity. This ramp maps the population in each band onto a
distinct momentum state. Measuring this momentum distribution quantifies band
populations in the moving lattice.

Transport along paths passing through the singular LBTP leads to interband
transitions that vary with the turning angle (Figs. 2, C and D). For
trajectories that enter the singularity and then reverse onto themselves
($\theta=0$), the population remains nearly entirely in the initial $n=1$
band. For trajectories that continue with constant tangent through the
singularity ($\theta=\pi$), the atoms undergo a near complete transition to
the upper $n=2$ band (seen also in Ref. (?)). Over the full range of $\theta$,
each population undergoes one cycle of oscillation.

The unit-cell wavefunction of the $n=1$ and $n=2$ Bloch states near the Dirac
point can be represented as a pseudo-spin-1/2 vector, with $s$-orbital Wannier
states at the lattice sites A and B representing the up- and down-spin basis
states. In this basis, the Bloch states are eigenstates of the Hamiltonian
$H_{\mathrm{LBTP}}=-\mathbf{B}(\mathbf{q})\cdot\boldsymbol{\sigma}$ where
$\mathbf{B}(\mathbf{q})$ is a pseudo-magnetic field that lies in the
transverse pseudo-spin plane and ${\boldsymbol{\sigma}}$ is the vector of
Pauli matrices. $\mathbf{B}(\mathbf{q})$ has a magnitude $B=\hbar
v_{\mathrm{g}}|\mathbf{q}-\mathbf{q}_{\mathbf{K}}|$ that varies linearly with
distance from the singularity, and as shown in Fig. 2B, it has an orientation
(in the proper gauge) that is radially outward from $\mathbf{K}$. Here,
$v_{\mathrm{g}}$ is the group velocity near the Dirac point. The $2\pi$
rotation of $\mathbf{B}(\mathbf{q})$ about the Dirac point is responsible for
the $\pi$-valued Berry phase of trajectories that encircle the Dirac point
(?).

Figure 2: Non-Abelian state rotations around a Dirac point. (A) Enlarged view
of Brillouin zone map. Atoms are loaded into the state
$n=1,\mathbf{q}=\mathbf{\Gamma}$ in the lattice, and transported (trajectory
marked by red arrows) at constant accelerations from
$\mathbf{\Gamma}\rightarrow\mathbf{K}$ and then from $\mathbf{K}$ to a final
point lying on a circle of diameter $0.8\norm{\mathbf{q_{K}}}$ centered at the
Dirac point. In this particular experiment, the atoms evolve between the $n=1$
and $n=2$ bands. Here, the color of the Brillouin zones does not indicate the
state of the atomic wavepacket along the trajectory. Rather, the color scheme
is used to interpret the band index of atoms after a band mapping measurement.
(B) Bloch states near the Dirac point are described as pseudo-spin-1/2 states
in an pseudo-magnetic field (black arrows) that points radially outward from,
and wraps once around, the Dirac point. (C) Band mapping images at the final
quasi-momentum, with overlain Brillouin zone maps, show the band population
vary with turning angle $\theta$. A third, short-length and adiabatic
translation step (not shown) ensures band mapping does not occur near a
Brillouin zone boundary so that there is no ambiguity in the band index of
atoms (see Section 3 of [SM]). The spatial scale is indicated by the black
scale bar of length 0.1 mm. (D) Fractional band populations
$N_{n}/N_{\mathrm{total}}$ ($n=1$: green, $n=2$: blue, sum of other bands:
gray) vs. $\theta$. Means and standard mean errors determined from 7 repeated
measurements. The green and blue dashed lines show a
$\mathrm{cos}^{2}(\theta/2)$ and $\mathrm{sin}^{2}(\theta/2)$ dependence,
respectively.

This pseudo-spin model explains our observations. The atomic pseudo-spin
entering the Dirac point along a ray experiences a pseudo-magnetic field whose
orientation $\mathbf{n}$ remains constant and whose magnitude smoothly tunes
to zero. Under this field, the initial-state pseudo-spin remains aligned along
$\mathbf{n}$. Departing the Dirac point, the pseudo-spin experiences a
magnetic field along a new orientation $\mathbf{m}$, with
$\mathbf{n}\cdot\mathbf{m}=\cos\theta$, and a magnitude increasing linearly
with time. The pseudo-spin is thus placed in a superposition of eigenstates,
with population $\cos^{2}(\theta/2)$ in the $\mathbf{m}$-oriented pseudo-spin
eigenstate ($n=1$ band) and $\sin^{2}(\theta/2)$ in the orthogonal state
($n=2$ band). This simple prediction is in good agreement with our data (Fig.
2D), with residual differences accounted by numerical simulations (Fig. S6) of
the dynamics of non-interacting atoms over the finite duration of our
experimental stages [test]The maximum duration of our experiment was limited
by the observed decay of atoms from the excited-energy Bloch states populated
during lattice acceleration..

We find that passage through the Dirac point produces a phase-coherent
superposition of band states. Such coherence is demonstrated by allowing the
atoms to evolve at the final point of the trajectory for a variable time
before measuring populations in a basis different from the local energy
eigenbasis. Temporal oscillations in these measurements, with a frequency
matching the calculated gap between the $n=1$ and $n=2$ bands, demonstrate the
coherence of the atomic state following transport (Fig. S3).

The energy-time uncertainty relation places a bound on how finely the singular
point can be located by our method. We consider a trajectory where the
acceleration has magnitude $a$ near the singularity. The system spends a time
$\delta t\sim(\hbar/ma)\delta q$ within $\delta q$ of the singularity; the
energy gap has magnitude $\delta E\sim\hbar v_{\mathrm{g}}\delta q$ in that
vicinity. Setting $(\delta t)(\delta E)\sim\hbar$ establishes that the band
structure is effectively gapless within a quasi-momentum distance of $\delta
q=R\sim\sqrt{ma\hbar/v_{\mathrm{g}}}$ of the singularity. That is, the
nonholonomy generated by the singular point should be observed also for finite
time trajectories that pass within the effective radius, $R$, of the
singularity.

We measure $R$ by driving the atoms along a family of trajectories, shown in
Fig. 3A, that connect between the initial $\mathbf{\Gamma}$-point to a final
$\mathbf{\Gamma}$-point that is one reciprocal lattice vector away, and
performing band mapping measurements at the final point. These trajectories
cross the boundary between the first and second Brillouin zones at nine
equally-spaced points along the $\mathbf{K^{\prime}}-\mathbf{M}-\mathbf{K}$
line. As shown in Fig. 3B for various traversal times, $\tau$, we observe that
trajectories that pass directly through either Dirac point yield a band
population distribution that is independent of $\tau$ with $\sim 3/4$ of the
atoms transferring to the upper band. In contrast, for traversal times that
are longer and for paths that veer farther from the Dirac points, the
transition between bands is increasingly suppressed, demonstrating that $R$
decreases with decreasing $a$ (increasing $\tau$). At first glance, there is
an analogy between the experiment of Fig. 3 and Majorana losses. The study
presented in Fig. 3 is quite close to Majorana’s description (?) and to the
picture of the “Majorana hole” that appears in a spherical quadrupole magnetic
trap. However, a difference between our experiment and magnetic traps is that
in traps the spin flips that occur because of transport through the state
geometry of the system lead to loss, whereas, in our experiment they lead to
transitions between trapped bands of the lattice.

Figure 3: Effective size of a Dirac singularity. (A) Illustration of quasi-
momentum trajectories over which the atoms are transported for the
measurements in B. Each trajectory connects two different
$\mathbf{\Gamma}$-points, but traverses one of nine equally spaced points
along the $\mathbf{K^{\prime}}-\mathbf{M}-\mathbf{K}$ line. (B) Fractional
$n=1$ band population plotted, for different trajectory traversal times
$\tau$, against the trajectory midpoints from (A). Means and standard mean
errors are generated from 3-5 repeated measurements. (C) An effective radius
$R$ (plotted normalized by $\hbar k=h/\lambda$) is defined for each $\tau$ by
the distance along the $\mathbf{K^{\prime}}-\mathbf{M}-\mathbf{K}$ line for
which the threshold $N_{1}/N_{\mathrm{total}}=0.5$ is fulfilled. $R$
diminishes with larger $\tau$.

The singularity at an LBTP can be characterized by two different experimental
methods: either by Berry phase measurements along trajectories that encircle
the singularity (?), or, as shown here, through state rotations produced along
trajectories that pass through the singularity. These two methods are related,
but nonequivalent. Berry phase measurements measure the integrated Berry flux,
which is determined from a diagonal element of the Berry connection matrix
${\bf A}^{nn}_{{\bf q}}$. A $\pi$-valued flux is found pinned to the singular
point. In contrast, the non-Abelian state rotations detected in our method
derive directly from the off-diagonal elements ${\bf A}^{nm}_{{\bf q}}$ with
$n=1$ and $m=2$ being the two crossing bands [SM]See supplementary materials
on Science Online. Further, different from Berry phase measurements, our
method can be regarded as measuring the Hilbert-Schmidt quantum distance
$d^{2}({\bf q},{\bf q}^{\prime})=1-|\langle u_{{\bf q}^{\prime}}^{1}|u_{{\bf
q}}^{1}\rangle|^{2}$ (?, ?) with $\mathbf{q}$ identified as a point along the
input path into – and $\mathbf{q}^{\prime}$ as a point along the exit path
from – the singularity. The oscillation of the $n=1$ band population as a
function of the turning angle reveals the quantum distance to undergo one
complete oscillation between zero and unity on a contour encircling the LBTP.

The distinction between these two methods is dramatic in the case of a QBTP.
Like the LBTP, a singular QBTP also carries concentrated Berry flux that is
restricted, assuming time-reversal and $C_{6}$ symmetry, to be 0 or $\pm 2\pi$
(?). However, these values of Berry phase are undetectable via interference
measurements. In contrast, our method uncovers the characteristic nonholonomy
of the singular QBTP and the concomitant modulation of the quantum distance
around the singular point.

For a singular QBTP, the geometric structure of Bloch states of the two
intersecting bands at the vicinity of the singularity can again be described
as those of a pseudo-spin in a pseudo-magnetic field. Different from the LBTP,
here, the pseudo-magnetic field, lying in the transverse pseudo-spin plane,
has a magnitude that increases quadratically with distance to the singularity,
and has an orientation that wraps by an angle of $4\pi$ along a path
encircling the singularity; see Fig. 4B.

Figure 4: Non-Abelian state rotations around a QBTP. (A) Atoms are prepared in
the $n=3$ band at $\mathbf{Q}$, transported (along red arrows) to the QBTP at
$\boldsymbol{\Gamma}$, and then transported to a final quasi-momentum for band
mapping. In this particular experiment, the atoms evolve between the $n=3$ and
$n=4$ bands. Here, the color of the Brillouin zones does not indicate the
state of the atomic wavepacket along the trajectory. Rather, the color scheme
is used to interpret the band index of atoms after a band mapping measurement.
(B) The pseudo-magnetic field (black arrows) describing the Bloch state
geometry wraps twice in orientation for one revolution around the QBTP. (C) A
plot of normalized band population as a function of $\theta$, collected by
analyzing band mapping data (Fig. S2). Red circles: $n=3$ band; purple
circles: $n=4$ band; gray circles: bands with $n\neq 3,4$. Means and standard
mean errors are determined from 7 repeated measurements. Dashed red (purple)
lines are predictions based on a simple pseudo-spin model for $n=3$ ($n=4$)
populations. Our numerical simulations suggest that the data does not reach
unity oscillation amplitude owing to nonadiabaticity in the band mapping
procedure (Fig. S6).

We probe this geometric structure at the QBTP that occurs at $\mathbf{\Gamma}$
between the $n=3$ and $n=4$ bands of the honeycomb lattice. For this, we first
load the Bose-Einstein condensate into the $n=3$ band of the lattice by
“inverse band mapping.” In previous work (?), such loading into excited bands
was realized with _moving_ atoms in a _static_ lattice; here, we realize
similar state-preparation with _static_ atoms loaded into a _moving_ lattice.
Specifically, we gradually increase the depth of a honeycomb lattice moving
with velocity $\mathbf{v}_{\mathbf{Q}}$ that, as shown in Fig. 4A, is located
within the third Brillouin zone in the extended zone scheme. Thence, we
accelerate atoms (in the lattice frame) at constant acceleration along the
path $\mathbf{Q}\rightarrow\mathbf{\Gamma}$, into the QBTP, and then, turning
by an angle $\theta$, accelerate the atoms at a different constant
acceleration out of the QBTP and to the edge of the Brillouin zone. Final
points are chosen as ones where populations in the $n=3$ and $n=4$ bands are
easily distinguished by band mapping.

We observe a nonholonomy at the QBTP that is distinct from that observed at
the LBTP. Specifically, we observe two cycles of oscillation in the final band
populations over the interval $\theta\in[0,2\pi]$. This behavior is explained
well by the pseudo-spin representation of the singular QBTP. Different from
the LBTP, here, pseudo-magnetic-field orientations along the incoming
($\mathbf{n}$) and outgoing ($\mathbf{m}$) paths are related as
$\mathbf{n}\cdot\mathbf{m}=\cos 2\theta$. The nonholonomy of the QBTP now
produces populations of $\cos^{2}(\theta)$ and $\sin^{2}(\theta)$ in the $n=3$
(initial) and $n=4$ bands, respectively.

The amplitude of the observed oscillation is lower than suggested by this
simple theory; again, we ascribe this difference to dynamical effects of our
finite-duration acceleration and band-mapping stages. Nevertheless, the
periodicity of the oscillations, combined with the known time-reversal and
$C_{6}$ symmetry of our lattice, unambiguously determines the topological
winding number around the QBTP to be well defined and equal to 2 [SM].

In conclusion, we have demonstrated transport of a quantum system through
singular band-touching points with different topological winding numbers. We
observe non-Abelian, coherent state rotation between bands. The dependence of
this rotation on the relative orientation of path tangents entering and
exiting the singular point unambiguously measures the winding number.

Our method of probing band structure could be applied to gain insight on other
band structure singularities and on interaction effects. It would be
interesting to study higher-order singular band-touching points, i.e., between
more than two bands. We deliberately minimize interaction effects in these
experiments, but in future work it will be interesting to observe potential
interaction-induced instabilities of Dirac points, QBTPs or other band
touching points. The path-dependent non-Abelian nonholonomy observed here may
also pertain to chemical systems, where potential energy surfaces are endowed
similarly with conical intersections (?, ?, ?, ?), suggesting a potential
route for quantum state control in optically driven molecules.

## References

  * 1. D. Xiao, M.-C. Chang, Q. Niu, Reviews of Modern Physics 82, 1959 (2010).
  * 2. A. Bansil, H. Lin, T. Das, Rev. Mod. Phys. 88, 021004 (2016).
  * 3. D. R. Yarkony, The Journal of Physical Chemistry A 105, 6277 (2001).
  * 4. G. A. Worth, L. S. Cederbaum, Annual Review of Physical Chemistry 55, 127 (2004).
  * 5. W. Domcke, D. R. Yarkony, Annual Review of Physical Chemistry 63, 325 (2012).
  * 6. M. S. Schuurman, A. Stolow, Annual Review of Physical Chemistry 69, 427 (2018).
  * 7. B. Simon, Phys. Rev. Lett. 51, 2167 (1983).
  * 8. M. V. Berry, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 392, 45 (1984).
  * 9. F. Wilczek, A. Zee, Phys. Rev. Lett. 52, 2111 (1984).
  * 10. A. Bohm, A. Mostafazadeh, H. Koizumi, Q. Niu, J. Zwanziger, The Geometric Phase in Quantum Systems (Springer Berlin Heidelberg, 2003).
  * 11. M. I. Katsnelson, K. S. Novoselov, A. K. Geim, Nature Physics 2, 620 (2006).
  * 12. Y. Zhang, Y.-W. Tan, H. L. Stormer, P. Kim, Nature 438, 201 (2005).
  * 13. L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, T. Esslinger, Nature 483, 302 (2012).
  * 14. L. Duca, et al., Science 347, 288 (2015).
  * 15. M. Tarnowski, et al., Phys. Rev. Lett. 118, 240403 (2017).
  * 16. T. Li, et al., Science 352, 1094 (2016).
  * 17. N. Fläschner, et al., Science 352, 1091 (2016).
  * 18. J.-W. Rhim, K. Kim, B.-J. Yang, Nature 584, 59 (2020).
  * 19. K. Sun, E. Fradkin, Phys. Rev. B 78, 245122 (2008).
  * 20. K. Sun, H. Yao, E. Fradkin, S. A. Kivelson, Physical Review Letters 103, 046811 (2009).
  * 21. J. M. Murray, O. Vafek, Physical Review B 89, 201110 (2014).
  * 22. J. Shah, S. Mukerjee, Physical Review B 103, 195118 (2021).
  * 23. E. McCann, V. I. Fal’ko, Physical Review Letters 96, 086805 (2006).
  * 24. E. McCann, M. Koshino, Reports on Progress in Physics 76, 056503 (2013).
  * 25. S. Pujari, T. C. Lang, G. Murthy, R. K. Kaul, Physical Review Letters 117, 086404 (2016).
  * 26. K. Hejazi, C. Liu, L. Balents, Physical Review B 100, 035115 (2019).
  * 27. K. Hejazi, C. Liu, H. Shapourian, X. Chen, L. Balents, Physical Review B 99, 035111 (2019).
  * 28. G. Jotzu, et al., Nature 515, 237 (2014).
  * 29. See supplementary materials on Science Online.
  * 30. The maximum duration of our experiment was limited by the observed decay of atoms from the excited-energy Bloch states populated during lattice acceleration.
  * 31. E. Majorana, Nuovo Cimento 9, 43 (1932).
  * 32. J. P. Provost, G. Vallee, Communications in Mathematical Physics 76, 289 (1980).
  * 33. Y. Hwang, J. Jung, J.-W. Rhim, B.-J. Yang, Physical Review B 103, L241102 (2021).
  * 34. T.-H. Leung, et al., Phys. Rev. Lett. 125, 133001 (2020).
  * 35. C. D. Brown, S.-W. Chang, M. N. Schwarz, T.-H. Leung, V. Kozii, A. Avdoshkin, J. E. Moore, D. Stamper-Kurn v1.1.0, Zenodo (2022); https://doi.org/10.5281/zenodo.6788172.
  * 36. A. Browaeys, et al., Phys. Rev. A 72, (2005).
  * 37. G.-F. Zhang, Y. Li, C. Wu, Physical Review B 90, 075114 (2014).
  * 38. C.-K. Chiu, A. P. Schnyder, Physical Review B 90, 205136 (2014).
  * 39. C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, S. Ryu, Reviews of Modern Physics 88, 035005 (2016).
  * 40. G. Montambaux, L.-K. Lim, J.-N. Fuchs, F. Piéchon, Physical Review Letters 121, 256402 (2018).
  * 41. P. Soltan-Panahi, et al., Nature Physics 7, 434 (2011).

## Acknowledgments

We thank E. Altman, M. Zaletel, N. Read and J. Harris for early insights into
this work. Funding: We acknowledge support from the NSF QLCI program through
grant number OMA-2016245 and also NSF grant PHY-1806362, and from the ARO
through the MURI program (grant number W911NF-17-1-0323). C. D. B.
acknowledges support from the National Academies of Science, Engineering, and
Medicine Ford Postdoctoral Fellowship program. V. K. and J. E. M. were
supported by the Quantum Materials program at LBNL, funded by the U.S.
Department of Energy under contract number DE-AC02-05CH11231. A. A. and J. E.
M. acknowledge support from the NSF under grant number DMR-1918065 and from a
Kavli ENSI fellowship. J. E. M. acknowledges support from a Simons
Investigatorship. Author contributions: All authors contributed substantially
to the work presented in this manuscript. C.D.B., S-W. C., M.N.S. and T-H. L.
acquired the data and maintained the experimental apparatus. C.D.B., S-W. C.,
and M. N. S. analyzed the data. C.D.B. prepared the manuscript. M. N. S., A.
A. and V. K. performed the numerical calculations. A. A. and V. K. performed
supporting theory calculations. J.E.M. and D. S-K. supervised the study. All
authors worked on the interpretation of the data and contributed to the final
manuscript. Competing interests: The authors declare no competing interests.
Data and materials availability: The data and the codes used for both data
analysis and theory curve generation are available at Zenodo C. D. Brown,
S.-W. Chang, M. N. Schwarz, T.-H. Leung, V. Kozii, A. Avdoshkin, J. E. Moore,
D. Stamper-Kurn v1.1.0, Zenodo (2022); https://doi.org/10.5281/zenodo.6788172.

## List of Supplementary Materials

Materials and Methods
Supplementary Text
Figs S1-S7
References (?, ?, ?, ?, ?, ?, ?)

## Materials and Methods

We prepare Bose-Einstein condensates of $1-2\times 10^{4}$ 87Rb atoms in the
the hyperfine state $\ket{F=1,m_{F}=-1}$ in an optical dipole trap with trap
frequencies $\omega_{x,y,z}\approx 2\pi\times\\{23,41,46\\}$ Hz. The chemical
potential of the harmonically trapped gas is roughly $0.4\,\mbox{kHz}\times
h$. We load a condensate into a honeycomb lattice with lattice depth $\approx
5E_{\mathrm{rec}}-10E_{\mathrm{rec}}$, where the photon recoil energy
$E_{\mathrm{rec}}/h=1/2\pi(\hbar/2m)*(2\pi/\lambda)^{2}=$ 2 kHz. The lattice
is formed by the mutual interference of three $\lambda=1,064$ nm wavelength
laser beams propagating in the horizontal plane, intersecting at equal angles
of $120^{\circ}$, and polarized in-plane. In ref. (?) the authors describe an
optical honeycomb lattice in which atoms experience a vector a.c. Stark shift
that varies across the two sites of the unit cell, causing an A-B site
asymmetry in potential depth. This Stark shift is also present in our system,
causing a small gap calculated to be 33 Hz at the K-point. Yet, we do not
observe, nor would we expect to observe such a small gap having a noticeable
effect on our transport experiment. In Fig. 3 of the main text, we show that
the effects of the Dirac point are seen even for trajectories that do not hit
the Dirac point, and that pass entirely through gapped regions of the band
structure, so long as the trajectory is traversed quickly with respect to the
residual gaps near the Dirac point. We expect that trajectories that probe the
Dirac point on timescales that are short compared to the inverse gap time
produced by the vector a.c. Stark shift will not be affected by that gap.

We operate our lattice in a “lattice of tubes” configuration, with only weak
confinement in the vertical direction. Our calculations suggest that, with the
lattice turned on, there are approximately 90 atoms per lattice site at the
center of the trap. One lattice beam has a fixed angular frequency $\omega$,
while the other two beams have time-dependent angular frequencies
$\omega_{1}(t)=\omega+\delta\omega_{1}(t)$ and
$\omega_{2}(t)=\omega+\delta\omega_{2}(t)$. Under a closed feedback loop that
dynamically varies the phase difference between the intersecting beams that
form the lattice, time-dependent relative detuning between lattice beams
allows us to accelerate the lattice potential reproducibly in two dimensions.
The lattice acceleration drives the atoms through a two-dimensional trajectory
of quasi-momentum in the lattice frame. After transporting the atoms along
this trajectory, we adiabatically ramp down the lattice potential to map Bloch
states onto a plane wave basis (band mapping) and then measure the atomic
velocity distribution using velocity-space focusing and resonant imaging. This
distribution reveals the lattice band populations prior to band mapping.

## Supplementary Text

In this supplementary material, we provide additional details of our
experiment and the theoretical framework used to understand our measurements.
In Section 1, we give a detailed description of the working principle of our
lattice translation setup. The time sequence that controls lattice
translations is then discussed in Section 2. In Section 3, we provide evidence
that the system undergoes a coherent evolution after passing a singularity. In
Section 4, we provide a complementary way to understand the results presented
in Fig. 3 of the main text. In Section 5, we characterize the effective size
of the QBTP between the third and fourth bands, with experiments that are
analogous to those carried out in Fig. 3 of the main text. In Section 6, we
provide a method for numerically simulating the time evolution of our system,
which is useful for choosing a set of parameters for experiments. Finally, to
shed light on the connection between our experimental results and the
topological properties of the system, we provide two-band models in Section 7
for both the LBTP and the QBTP that reproduces the non-Abelian state rotations
that we measure.

###### Contents

  1. 1 Arbitrary lattice translation
  2. 2 Lattice translation sequence
  3. 3 Coherent evolution after population transfer at a singularity
  4. 4 Complementary representation of data in Fig. 3 of the main text
  5. 5 Momentum space extent of the QBTP
  6. 6 Simulations via full Hamiltonian construction
  7. 7 Population dynamics and topology of the band touching points within the two-band models
    1. 7.1 Linear band touching
    2. 7.2 Quadratic band touching

##  1 Arbitrary lattice translation

The optical honeycomb lattice potential is created by interfering three laser
beams, with wavelength $\lambda=1,064$ nm, that intersect at relative angles
of $120^{\circ}$. Written explicitly, the potential is given by

$\displaystyle
V=\left|\sqrt{\frac{2V_{\mathrm{lat}}}{9}}\sum_{i=1,2,3}\mathbf{\epsilon}_{i}e^{\mathrm{i}\mathbf{k}_{i}\cdot\mathbf{r}}\right|^{2}=\frac{2}{9}V_{\mathrm{lat}}\left(3-\sum_{i=1,2,3}\cos{(\mathbf{G}_{i}\cdot\mathbf{r})}\right),$
(S1)

where $\mathbf{k}_{n}=k(\cos{\theta_{i}},\sin{\theta_{i}},0)$ 111Here we use a
Cartesian basis with $\mathbf{x}$ and $\mathbf{y}$ being unit vectors spanning
the horizontal lattice plane, and $\mathbf{z}$ oriented vertically. are the
wavevectors of lattice beams with
$k=2\pi/\lambda,\,\theta_{i}=(-1+\frac{2i}{3})\pi$,
$\mathbf{\epsilon}_{i}=(\sin{\theta_{i}},-\cos{\theta_{i}},0)$ are the (in-
plane) polarization vectors, $V_{\mathrm{lat}}$ is the depth of the lattice,
and $\mathbf{G}_{i}=\sum_{j,k}\epsilon_{ijk}\mathbf{k}_{j}$ are the reciprocal
lattice vectors ($\epsilon_{ijk}$ is the antisymmetric Levi-Civita tensor).
Here, the script index $i$ is not to be confused with the imaginary number
$\mathrm{i}$. When one of the beams with angular frequency
$\omega_{1}=\omega+\delta\omega_{1}$, is detuned by $\delta\omega_{1}$
relative to the beam with angular frequency $\omega$, the lattice will
translate at a velocity $\mathbf{v}_{1}=a\Delta\omega_{1}\hat{k}_{1}$ in the
lab frame, where $a=\frac{2\lambda}{3}$ is the lattice constant. Here, we
ignore the small change in the wavevector $\mathbf{k}_{1}$ due to the
detuning. Similarly, detuning beam 2 by $\Delta\omega_{2}$ translates the
lattice at a velocity of $\mathbf{v}_{2}=a\Delta\omega_{2}\hat{k}_{2}$.
Therefore, by choosing an appropriate linear combination of $\mathbf{v}_{1}$
and $\mathbf{v}_{2}$, we can translate the lattice at any velocity in the
lattice plane. Acceleration along a linear trajectory in quasi-momentum space
corresponds to varying the detunings of the two lattice beams linearly in
time.

In the experimental setup, each lattice beam is controlled by an acousto-optic
modulator (AOM), which we use to control the detuning between lattice beams.
We change the drive frequency of AOMs by tuning the set voltage of the
voltage-controlled oscillators (VCOs) that control the AOMs, as shown in Fig.
S1B. For each line segment in momentum space, we change the detuning linearly,
which creates a uniform acceleration of the lattice. In the reference frame
where the lattice is stationary, the atoms are thus accelerated in the
opposite direction.

Figure S1: Experimental sequence. (A) Typical ramp shape of lattice beam
intensity. The same ramp function is shared for all three lattice beams. (B)
Typical ramp shape of lattice beam detuning. Generally, the detunings between
two beam pairs are different, depending on the trajectory in momentum space.
The shaded bars on top and bottom of both A and B indicate the time intervals
(in units of milliseconds) used for each stage of the experiment: top (gray) =
slow load, bottom (red) = fast load.

##  2 Lattice translation sequence

The ramp shapes of beam intensities are given in Fig. S1A. After creating a
Bose-Einstein condensate, we first adiabatically ramp up the lattice
intensity, such that all atoms are loaded into the ground band at
$\mathbf{\Gamma}$; see “slow load” in Fig. S1. This also allows the atoms to
follow the minimum of the combined trap potential created by the optical
dipole trap (ODT) beams and the lattice beams, which varies as the lattice
intensity is ramped up, due to possible misalignment between ODT beams and
lattice beams.

For the QBTP experiment we prepare atoms in the $n=3$ band, which can be done
by rapidly loading atoms into a running lattice. In reality, to make sure that
the atoms are at the minimum of the trap potential at the start of the
experiment, we first load atoms into the ground band with the slow loading
procedure. Next, we perform band mapping to recover a stationary BEC, before
rapidly ramping up the lattice beam intensities again but with a running
lattice. This is shown in Fig. S1 as “fast load.” With this technique, we can
load into the $n=3$ band with $>92\%$ fidelity on average.

Band population detection (or “band mapping”) is achieved by ramping down the
lattice beams slowly enough, such that atoms in different bands are
adiabatically mapped to plane waves with momenta separated by reciprocal
lattice vectors. In the resulting absorption images, the positions of atom
peaks then correspond to their momentum in the lab frame. In the main text, we
show our band mapping measurement results in the lattice frame, which is
obtained by translating the lab frame image by a distance determined by the
final quasi-momentum in a lattice acceleration sequence. The positions of atom
peaks in the images now correspond to their momenta in the lattice frame, and
the association between band label $n$ and these peaks can then be determined
by overlapping the Brillouin zone of the lattice with the image. Fig. S2
provides an example of images obtained from band mapping measurements.

In the measurement of Fig. 2C of the main text, if the angle of the outgoing
trajectory would lie close to a Brillouin zone edge, we accelerate the atoms
further before band mapping to a point far away from any Brillouin zone edges.
This further acceleration is slow and adiabatic, such that the measured
population in each band is not affected by the additional translation.

Figure S2: Band mapping data. Shown in the figure are individual measurements
from our characterization of the honeycomb-lattice QBTP, as summarized in Fig.
4 of the main text. Images $\mathrm{i}$ through $\mathrm{iv}$ (in A and C)
show the results of band mapping after trajectories with four different
turning angles $\theta$, at the final quasi-momenta marked by black points on
the Brillouin zone map in B. (A) Band mapping images obtained in the lab
frame. These are absorption images of the atomic spatial distribution taken
after (1) adiabatically ramping off the lattice beams (and also optical traps)
at the final lattice velocity, and (2) allowing the atoms to expand in a weak
harmonic magnetic confinement so as to map velocity onto position. The spatial
scale of the images in A and C is indicated by the black scale bar of length
0.1 mm. The scale of the images in terms of velocity distributions is
indicated by the overlain Brillouin zone maps. The color bar applies to both A
and C. (C) The band mapping images in the lattice frame. Each of them are
obtained by translating the lab-frame image in A by the corresponding quasi-
momentum (one of the red arrows) in B.

##  3 Coherent evolution after population transfer at a singularity

Here, we show that after the population transfer at $\mathbf{K}$, the atoms in
bands $n=1$ and $n=2$ still interfere coherently. The trajectory used for the
coherence test is shown in Fig. S3A. After creating an even superposition in
the lowest two bands by turning at $\mathbf{K}$, the atoms are accelerated to
the center of the second Brillouin zone, at point $\mathbf{Q}$ in momentum
space, and held there for a varying amount of time $t$. During the hold, the
state acquires a relative phase $\varphi=t\Delta\omega_{\mathbf{Q}}$ between
the $\ket{n=1}$ and $\ket{n=2}$ components, where
$\hbar\Delta\omega_{\mathbf{Q}}$ is the band gap between the lowest two bands
at $\mathbf{Q}$. Then, we accelerate the atoms quickly back to
$\mathbf{\Gamma}$, with an acceleration large enough such that the state is
projected onto the basis states formed by the Bloch states at
$\mathbf{\Gamma}$. Finally, we perform a band mapping measurement to project
the phase difference into a population difference between bands $n=1$ and
$n=2$. In Fig. S3B, we indeed see an oscillation in the measured population
that is consistent (with the empirically added decay) with what is expected
from theory.

Figure S3: Detecting coherence of superposition produced by transport through
the LBTP singularity. (A), Atoms are transported on a trajectory with three
constant accelerations: $\boldsymbol{\Gamma}\rightarrow\mathbf{K}$,
$\mathbf{K}\rightarrow\mathbf{Q}$, and
$\mathbf{Q}\rightarrow\boldsymbol{\Gamma}$. The first two acceleration steps
are done adiabatically, in 0.3 ms each. Atoms are then held at $\mathbf{Q}$
for a variable hold time. The last acceleration step is done quickly in 0.2
ms. This faster acceleration causes the atomic state (a superposition in the
$n=1$ and $n=2$ bands at $\mathbf{Q}$) to be projected onto the band
eigenbasis at $\boldsymbol{\Gamma}$ prior to band mapping. (B) $n=1$ band
population seen from band mapping at $\boldsymbol{\Gamma}$. The dynamical
phase acquired by the superposition state $\mathbf{Q}$ leads to temporal
oscillation in the band-mapping $n=1$ population. We fit the data to a
sinusoid with a decaying exponential envelope (red curve), which finds that
the frequency is $1.47\pm 0.53$ kHz, and the decay constant $319\pm 205\
\mathrm{\mu}$s. This frequency is in agreement with the band gap between $n=1$
and $n=2$ bands at $\mathbf{Q}$, which is calculated as $1.77$ kHz at a
lattice depth of $V_{\mathrm{lat}}=h\times 24$ kHz $=12E_{\mathrm{rec}}$ for
non-interacting atoms. The black curve is the simulation result obtained by a
numerical simulation of the entire transport, hold time, and projective
measurement process, in the non-interacting limit, as described in Sec. 6. The
simulation does not account for decay: To make a more direct comparison
between theory and experiment, we apply exponential damping to the simulation
result using the decay constant found from our fit to experimental data.

##  4 Complementary representation of data in Fig. 3 of the main text

In Fig. 3 of the main text, we show a measurement of the effective size of a
Dirac point singularity by plotting the population in the lowest band against
the midpoint of various trajectories (the trajectories are shown both in Fig.
3 of the main text, and in Fig. S4A). Instead, here we provide a different
visualization of the data. In Fig. S4B we plot the $n=1$ band population
normalized by the total atom number against the acceleration time used to
traverse the trajectories in Fig. S4A. A projection of the initial state
$\ket{\psi^{1}_{\mathbf{\Gamma}}}$ onto the $n=1$ Bloch state at the final
$\mathbf{\Gamma}$-point in the two-band model gives 0.25, which is the
expected normalized population in the $n=1$ band under the sudden
approximation. For trajectories far away from singularities, the adiabatic
theorem predicts that the $n=1$ band population should increase with longer
acceleration times, which we indeed observe in our experiments. On the other
hand, when we cross $\mathbf{K}$ in the trajectory, the band gap between $n=1$
and $n=2$ bands vanishes, and the adiabatic theorem no longer applies.

Figure S4: Complementary representation of the same data set of Fig. 3 in the
main text. (A) Trajectories taken in the experiment copied here for reference.
(B) normalized $n=1$ band population as a function of acceleration time
$\tau$. Only data for trajectories that go through points between $\mathbf{K}$
and $\mathbf{M}$ are shown in the plot, for simplicity. Here, the solid lines
are the simulation results obtained by numerical simulations of the entire
transport, hold time, and projective measurement process, in the non-
interacting limit, as described in Sec. 6. For the shortest acceleration time
used, we expect the $n=1$ band population to deviate from 0.25 (dashed line),
predicted by the two-level model in Sec. 7, because of excitation to higher
bands.

##  5 Momentum space extent of the QBTP

To characterize the singular behavior of the QBTP, we perform a measurement
similar to the one shown in Fig. 3 of the main text. As seen in Fig. S5B, for
the same acceleration time $\tau$, more atoms stay in the third band for the
trajectories that are farther away from $\mathbf{\Gamma}$, since the condition
for adiabatic evolution is better satisfied. At very short acceleration times,
transitions to higher bands are significant, which is also captured by
simulations. At longer acceleration times ($\tau>1.4$ ms), the actual quasi-
momentum of atoms in the lattice deviates from the set values due to
acceleration in the lab frame, which appears in the images as an overall shift
of peak positions after band mapping (Fig. S5C). There are two immediate
consequences: first, the actual trajectory is different from expected; second,
band mapping becomes unreliable if the actual final quasi-momentum is close to
the boundary between the third and fourth Brillouin zones. Therefore, we only
plot the measurement results up to $\tau=2.8$ ms, after which we have
difficulties assigning band indices to different wave packets in the images.
We attribute this acceleration to the fact that atoms are dragged along the
direction of lattice translation during the sequence (plus corrections due to
group velocity of atoms in the lattice) (?). The dragging effect grows with
the duration of the experiment. The QBTP experimental sequence is
significantly longer than the LBTP experiment (there are more steps in the
procedure and the band mapping time is necessarily longer), so the dragging
effect is more pronounced. Once the atoms are dragged away from the potential
minimum of the optical trap, this confining potential applies an acceleration
to the atoms. Finally, we indeed observe that, like the LBTP, the region in
which band population rotations occur is confined at the QBTP.

Figure S5: Momentum-space extent of the QBTP. (A) Trajectories chosen for the
experiments. (B) The population in each band as a function of acceleration
time $\tau$. Here, $\tau$ is the total time taken to accelerate the atoms from
the initial state to the final state. For this choice of band mapping point,
only the third and fourth bands are uniquely distinguished. Error bars are one
standard error. (C), Acceleration in lab frame. As described in Sec. 2, the
positions of atom peaks in images correspond to their lab frame velocities.
Therefore, by tracing the position of one of the peaks in band mapping images
in the lab frame, we can measure the change in lab-frame momentum that atoms
obtain during the experiment. The figure is plotted with data from B-i, and
different colors are assigned for each acceleration time labeled next to the
data points. The radius of the dashed circle is $0.1\|\mathbf{q_{K}}\|$,
centered at the position corresponding to zero momentum in lab frame, which is
determined by repeated measurements of the atomic wavepacket’s position
without lattice translation in a separate set of experiments.

##  6 Simulations via full Hamiltonian construction

The dynamics of the system is governed by the time dependent Schrödinger-
equation:

$\mathrm{i}\hbar\partial_{t}\ket{\psi}=H(t)\ket{\psi}.$ (S2)

We divide the time sequence into small time steps with effectively constant
Hamiltonian $H(t)$ such that we can use the time-independent Schrödinger
equation. Bloch’s theorem allows us to write the Hamiltonian in the plane wave
basis
$\innerproduct{\mathbf{r}|\mathbf{G}_{(n,m)}}{\mathbf{r}|\mathbf{G}_{(n,m)}}=\exp{(\mathrm{i}\mathbf{G}_{(n,m)}\cdot\mathbf{r})}$,
where
$\mathbf{G}_{(n,m)}=n\,\mathbf{G}_{1}+m\,\mathbf{G}_{2}\;(n,m\in\mathbb{Z})$
is a reciprocal lattice vector of the $1,064\,$nm 2D optical honeycomb
lattice. The Schrödinger equation then takes the form

$\bigg{(}\frac{1}{2m}\big{(}\mathbf{p}+\hbar\mathbf{q}\big{)}^{2}+V(\mathbf{r})\bigg{)}u_{\mathbf{q}}^{n}(\mathbf{r})=E_{\mathbf{q}}^{n}u_{\mathbf{q}}^{n}(\mathbf{r}),$
(S3)

where $\mathbf{q}$ is the quasi-momentum vector, $n$ is the band index and
$u_{\mathbf{q}}^{n}(\mathbf{r})=\langle{\bf r}|u^{n}_{\bf q}\rangle$ is the
real space representation of the lattice periodic part $|u^{n}_{\bf q}\rangle$
of a Bloch state
$\ket{\psi_{\mathbf{q}}^{n}}=e^{\mathrm{i}\mathbf{q\cdot\hat{r}}}|u^{n}_{\bf
q}\rangle$. By expanding $V(\mathbf{r})$ and $u_{\mathbf{q}}^{n}(\mathbf{r})$
in the plane wave basis the Hamiltonian can be expressed explicitly as

$H_{(n,m),(n^{\prime},m^{\prime})}(\mathbf{q})=\bigg{(}\frac{2}{3}V_{0}+\frac{\hbar^{2}k^{2}}{2m}|\tilde{\mathbf{G}}_{(n,m)}+\tilde{\mathbf{q}}|^{2}\bigg{)}\delta_{n,n^{\prime}}\delta_{m,m^{\prime}}-\frac{V_{0}}{9}\delta_{n,n^{\prime}+j}\delta_{m,m^{\prime}+l}\;\;\text{with
}(j,l)\in\mathcal{J},$ (S4)

where $\mathcal{J}=\\{(-1,0),(0,-1),(1,0),(0,1),(-1,-1),(1,1)\\}$, the tilde
implies that the quantity is normalized by the lattice laser wavenumber
$k=2\pi/\lambda$ with $\lambda=1,064~{}\mathrm{nm}$ and $V_{0}$ as the optical
lattice depth. Here, $V_{0}<0$ because the lattice light is red-detuned with
respect to the principal atomic resonances of 87Rb.

We simulate the system dynamics by truncating the plane wave basis such that
it incorporates 81 different plane wave components and we divide the time
evolution into time steps of $\tau=5~{}\mu$s. We ensured that the simulation
captures all relevant dynamics of the system with these settings. The initial
state of the system $\ket{\psi(0)}$ is often taken to be the ground band
($n=1$) eigenstate. We utilize spectral decomposition to speed up our
computation. We apply $H(t)$ to the state such that
$\ket{\psi(t)}=T\prod_{t^{\prime}<t}e^{-\mathrm{i}H(t^{\prime})\tau/\hbar}\ket{\psi(0)}$,
where $T$ implies that the sum is properly time ordered. With this method we
can easily extract the populations of the state in each of the eigenstates
$|u^{n}_{\bf q}\rangle$ at each time step by first calculating the eigenstates
and subsequently the overlap with the state $\ket{\psi(t)}$, as shown in Fig.
S6A. In that manner we can ensure that the effects we observe originate in the
system dynamics and that we are able to quantify the effect of non-adiabatic
population transfer. We use this technique to minimize excitations due to
dynamics unrelated to the topological nature of the band structure
singularities.

Figure S6: Simulation by decomposition of dynamic system state into Bloch
eigenfunctions of $H(t)$. (A) Example of state evolution for
$\theta=255^{\circ}$ in Fig. 4 of the main text. The state
$\ket{\psi(0)}=\ket{\mathbf{\Gamma},1}$ evolves over the following steps: i)
lattice ramp down, ii) ramp to $\mathbf{q}=1.25\,\mathbf{K}$, iii) rapid
lattice intensity ramp up, iv) ramp to $\mathbf{q}=\mathbf{\Gamma}$, v) ramp
to $\mathbf{q}=\mathbf{q}_{255\text{\textdegree}}^{r=0.9}$, before vi) band
mapping. In v), the superscript labels the distance away from
$\mathbf{\Gamma}$, in units of $\mathbf{q_{K}}$. Due to nonadiabatic
evolution, bands other than $n=3,4$ are populated. (B) and (C) Simulation
results (squares with dashed lines) for the experiments carried out in Fig. 2
and Fig. 4, respectively. Experiment data (circles) are copied here for
reference.

The simulated population transfer results for the experiments in Fig. 2 and
Fig. 4 of the main text are plotted in Fig. S6B and Fig. S6C, with the
experiment data reproduced in the same plot. Although these simulations do not
involve any interaction effects, they nevertheless show no significant
discrepancies to the observed data.

##  7 Population dynamics and topology of the band touching points within the
two-band models

In this section, we show that all the main results presented in the main text
can be qualitatively understood within an effective two-band model. To do
that, we derive the equations which describe the band population of a particle
moving in a lattice in the presence of an externally applied force within a
two-band model. Relevant to our experiment, we only focus on the case of a
constant force. In our derivation, we closely follow that from the
supplementary materials from Ref. (?). Note that in the following sections $i$
is the imaginary number.

We start with the two-band lattice Hamiltonian of the form

$\hat{H}_{0}=\sum_{n=1,2}\sum_{{\bf q}}{\varepsilon}_{{\bf q}}^{n}|\psi_{{\bf
q}}^{n}\rangle\langle\psi_{{\bf q}}^{n}|,$ (S5)

where $|\psi_{{\bf q}}^{n}\rangle$ are Bloch eigenstates of the $n$-th band at
quasi-momentum ${\bf q}$ having energy ${\varepsilon}_{{\bf q}}^{n}$. The
Bloch states can be parameterized as $|\psi^{n}_{\bf q}\rangle=e^{i{\bf
q}\cdot\hat{\bf r}}|u^{n}_{\bf q}\rangle$, where $|u^{n}_{\bf q}\rangle$ has
periodicity of the lattice in real space and $\hat{\bf r}$ is the position
operator. The Bloch states satisfy the normalization conditions
$\langle\psi^{n}_{\bf q}|\psi^{m}_{{\bf q}^{\prime}}\rangle=\delta({\bf
q}-{\bf q}^{\prime})\delta_{nm}$ and $\langle u^{n}_{\bf q}|u^{m}_{\bf
q}\rangle=\delta_{nm}$.

In the presence of an external constant force ${\bf F}$ the Hamiltonian of the
system reads as $\hat{H}=\hat{H}_{0}-{\bf F}\cdot\hat{\bf r}$, leading to the
Schrödinger equation

$i\partial_{t}|\Psi(t)\rangle=(\hat{H}_{0}-{\bf F}\cdot\hat{\bf
r})|\Psi(t)\rangle,$ (S6)

where we use the units with $\hbar=1$.

We assume that at $t=0$ the particle has quasi-momentum ${\bf q}_{0}$,
implying that its wavefunction is given by

$|\Psi(0)\rangle=\sum_{n=1,2}\varphi_{n}(0)|\psi_{{\bf q}_{0}}^{n}\rangle,$
(S7)

and the band population at $t=0$ equals $|\varphi_{n}(0)|^{2}$. In the
presence of an external force, the quasi-momentum of the particle changes
according to

${\bf q}(t)={\bf q}_{0}+{\bf F}t,$ (S8)

meaning that the solution has form

$|\Psi(t)\rangle=\sum_{n=1,2}\varphi_{n}(t)|\psi_{{\bf q}(t)}^{n}\rangle.$
(S9)

The Schrödinger equation can then be rewritten as:

$i\partial_{t}\left(\begin{array}[]{c}\varphi_{1}\\\
\varphi_{2}\end{array}\right)=\left(\begin{array}[]{cc}{\varepsilon}^{1}_{{\bf
q}(t)}-\xi^{11}_{{\bf q}(t)}&-\xi^{12}_{{\bf q}(t)}\\\ -\xi^{21}_{{\bf
q}(t)}&{\varepsilon}^{2}_{{\bf q}(t)}-\xi^{22}_{{\bf
q}(t)}\end{array}\right)\left(\begin{array}[]{c}\varphi_{1}\\\
\varphi_{2}\end{array}\right),$ (S10)

where we have defined

$\xi^{nm}_{{\bf q}(t)}\equiv{\bf F}\cdot{\bf A}^{nm}_{{\bf q}(t)}=i\langle
u^{n}_{{\bf q}(t)}|\partial_{t}|u^{m}_{{\bf q}(t)}\rangle,$ (S11)

and

${\bf A}^{nm}_{{\bf q}}\equiv i\langle u^{n}_{{\bf q}}|\partial_{{\bf
q}}|u^{m}_{{\bf q}}\rangle$ (S12)

is the Berry connection. As we see from Eq. (S10), it is precisely the off-
diagonal components of the Berry connection that drives the population change
of the two bands.

###  7.1 Linear band touching

To study linear band touching, we apply the above analysis to the two lowest
bands of the typical band structure of the honeycomb lattice arising from the
$s$-orbitals. Under the tight-binding approximation with only nearest neighbor
hopping, the Hamiltonian of the (non-driven) system has the form

$H_{0}({\bf q})=\left(\begin{array}[]{cc}0&M({\bf q})\\\ M^{*}({\bf
q})&0\end{array}\right),$ (S13)

with

$M({\bf
q})=-t\left(2e^{iq_{y}a/2}\cos\frac{\sqrt{3}q_{x}a}{2}+e^{-iq_{y}a}\right),$
(S14)

where $-t$ is the hopping integral and $a$ is the lattice constant.

With this model Hamiltonian, we use Eq. (S10) to calculate the final band
population as the particles move along the trajectories shown in Fig. S4A (or
Fig. 3A of the main text), assuming that originally they were prepared in the
$n=1$ state. The result as a function of the trajectory midpoint and the
acceleration time $\tau$ is presented in Fig. S7 and is in good qualitative
agreement with Figs. 3B,C of the main text and Fig. S4B. The only difference
appears at short acceleration times, which is quite natural since at such
large acceleration the population of the higher bands becomes substantial and
must be taken into account.

Figure S7: Final band population as a function of the acceleration time $\tau$
and different trajectories shown in Fig. S4A calculated within a two-band
model, Eq. (S13). Left: Population of the $n=1$ band as a function of the
trajectory midpoint at different fixed $\tau$ and the size of the singularity
$R$ as a function of $\tau$. The black arrow indicates the increase of $\tau$.
The results are in good qualitative agreement with the measured data shown in
Figs. 3B,C of the main text. Right: The normalized $n=1$ band population
presented as a function of the acceleration time for different fixed
trajectories. Again, the result of the calculation is in good qualitative
agreement with the experimental data shown in Fig. S4B.

The Dirac singularities are located at the points
$\mathbf{q}_{\mathbf{K^{\prime}}}=\left(4\pi/3\sqrt{3}a,0\right)$ and
$\mathbf{q}_{\mathbf{K}}=\left(2\pi/3\sqrt{3}a,2\pi/3a\right)$ in the
Brillouin zone and those related by the reciprocal lattice vectors. Near the
$\mathbf{K}$-point, the Hamiltonian can be expanded to linear order in quasi-
momentum ${\bf k}\equiv{\bf q}-\mathbf{q}_{\mathbf{K}}$ as

$H_{0}({\bf k})\approx v_{g}\left(\begin{array}[]{cc}0&k_{+}\\\
k_{-}&0\end{array}\right),\qquad k_{\pm}=k_{x}\pm ik_{y},$ (S15)

where $v_{g}=3at/2$ is the group velocity at the Dirac point. The wave
functions of the $n=1$ and $n=2$ bands are given by

$|u_{\bf k}^{1}\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}-k_{+}/k\\\
1\end{array}\right),\qquad|u_{\bf
k}^{2}\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}k_{+}/k\\\
1\end{array}\right),$ (S16)

with $k_{\pm}\equiv k_{x}\pm ik_{y},$ $k\equiv|{\bf
k}|=\sqrt{k_{x}^{2}+k_{y}^{2}}$, and the components of the Berry connection
equal

${\bf A}_{{\bf k}}^{11}={\bf A}_{{\bf k}}^{22}=-{\bf A}_{{\bf k}}^{12}=-{\bf
A}_{{\bf
k}}^{21}=\frac{1}{2(k_{x}^{2}+k_{y}^{2})}\left(\begin{array}[]{c}k_{y}\\\
-k_{x}\end{array}\right).$ (S17)

We emphasize that while the specific form of the wave functions and Berry
connection depend on the gauge of the wave functions, the final answer for the
observables, e.g. band population, is gauge-independent.

When a particle moves exactly towards or away from the singularity in momentum
space, it does not experience any interband transitions away from the singular
point itself since ${\bf F}\cdot{\bf A}_{{\bf k}}^{mn}=0$. This is a general
feature of any homogeneous Hamiltonian satisfying the relation $H_{0}(s{\bf
k})=s^{\alpha}H_{0}({\bf k})$. Indeed, if choosing a proper gauge like the one
in Eq. (S16), the eigenstates of the Hamiltonian depend only on the direction
of the quasi-momentum, $\hat{\bf k}={\bf k}/|{\bf k}|$, but not on its
absolute value. Then, according to Eq. (S11), ${\bf F}\cdot{\bf A}_{{\bf
k}}^{mn}=0$ because $\hat{{\bf k}}(t)$ and, consequently, $|u_{{\bf
k}(t)}^{m}\rangle$ do not depend on $t$. This implies that all the transitions
happen exactly at the singularity. This conclusion is true provided we stay in
the region of the Brillouin zone where the Hamiltonian of the system can be
approximated by a homogeneous function (i.e., we must remain sufficiently
close to the singular point). Assuming that the particle moving along the
direction $\hat{\bf k}$ scatters off the singularity in the direction
$\hat{\bf k}^{\prime}$ at time $t_{0}$, the transition probability between the
two bands can then be found from the continuity condition
$|\Psi(t_{0}-0)\rangle=|\Psi(t_{0}+0)\rangle$, which translates into

$\displaystyle\varphi_{1}(t_{0}-0)|u_{\bf
k}^{1}\rangle+\varphi_{2}(t_{0}-0)|u_{\bf k}^{2}\rangle$ $\displaystyle=$
$\displaystyle\varphi_{1}(t_{0}+0)|u_{{\bf
k}^{\prime}}^{1}\rangle+\varphi_{2}(t_{0}+0)|u_{{\bf k}^{\prime}}^{2}\rangle.$
(S18)

If the particle was originally in the $n=1$ band, ${\varphi}_{2}(t_{0}-0)=0$,
then the transition probability is given by

$d^{2}({\bf k},{\bf k}^{\prime})=|{\varphi}_{2}(t_{0}+0)|^{2}=|\langle u_{{\bf
k}^{\prime}}^{2}|u_{{\bf k}}^{1}\rangle|^{2}=1-|\langle u_{{\bf
k}^{\prime}}^{1}|u_{{\bf k}}^{1}\rangle|^{2}.$ (S19)

This quantity is exactly the Hilbert-Schmidt quantum distance between the
states $|u_{{\bf k}^{\prime}}^{1}\rangle$ and $|u_{{\bf k}}^{1}\rangle$ (?,
?). For two infinitesimally close states ${\bf k}$ and ${\bf k}^{\prime}={\bf
k}+\mathrm{d}{\bf k}$, it reveals the gauge-invariant quantum metric tensor
$g_{ij}({\bf k})$:

$d^{2}({\bf k},{\bf k}+d{\bf k})=g_{ij}({\bf
k})\mathrm{d}k_{i}\mathrm{d}k_{j},$ (S20)

which is defined as

$g_{ij}({\bf k})=\frac{1}{2}[\langle\partial_{i}u_{\bf k}|\partial_{j}u_{\bf
k}\rangle-\langle\partial_{i}u_{\bf k}|u_{\bf k}\rangle\langle u_{\bf
k}|\partial_{j}u_{\bf k}\rangle+(i\leftrightarrow j)].$ (S21)

In the case of the isotropic Dirac dispersion, the quantum distance is
expressed through the angle $\theta$ between the momenta ${\bf k}$ and ${\bf
k}^{\prime}$:

$d^{2}_{\text{Dirac}}({\bf k},{\bf k}^{\prime})=\sin^{2}\frac{\theta}{2}.$
(S22)

This result is in perfect agreement with the pseudo-magnetic field picture
discussed in the main text and leads to the same prediction.

It is worth noting that such a clear and simple relation between the band
population change and the quantum distance (and Hilbert space geometry) is,
strictly speaking, only possible in the vicinity of the singular point, where
the Hamiltonian can be approximated by a homogeneous function. In the
experiment, on the other hand, the start and end points of the trajectory in
the Brillouin zone clearly lie outside the singular region. Nevertheless, all
the conclusions derived above remain valid provided the particle’s motion is
sufficiently slow. Indeed, outside the singular region, the system is
effectively gapped. Hence, according to the adiabatic theorem, the interband
transitions are rare if the drive is slow.

###  7.2 Quadratic band touching

The entire analysis for the quadratic band touching point is similar to that
of the previous section, yet the results show some important differences,
which we highlight and discuss in detail below.

The lowest-energy quadratic band touching is formed by the $n=3$ and $n=4$
bands at the $\mathbf{\Gamma}$-point. To derive the effective low-energy
theory for the QBTP, we follow Refs. (?, ?) and start with the tight-binding
model band structure that arises from the $p_{x}$ and $p_{y}$ orbitals. We
only include the $\sigma$-bonding and nearest-neighbor hopping for simplicity.
The Hamiltonian in this case takes the form

$H=t_{\sigma}\sum_{i=1}^{3}\sum_{\vec{r}\in
A}\left\\{p^{\dagger}_{i}(\vec{r})p_{i}(\vec{r}+a\hat{e}_{i})+\text{H.
c.}\right\\},$ (S23)

where $p_{i}=(p_{x}\hat{e}_{x}+p_{y}\hat{e}_{y})\cdot\hat{e}_{i}$ are the
projections of the $p$ orbitals parallel to the nearest-neighbor bond
directions $\hat{e}_{i}$, $p_{x,y}^{\dagger}$ represent the creation operators
of atoms on the corresponding orbitals, and $t_{\sigma}$ is the
$\sigma$-bonding strength. Vectors $\hat{e}_{i}$ are three unit vectors from
one $A$ site to its three neighboring $B$ sites and are given by
$\hat{e}_{1,2}=\pm\frac{\sqrt{3}}{2}\hat{e}_{x}+\frac{1}{2}\hat{e}_{y}$ and
$\hat{e}_{3}=-\hat{e}_{y}$.

The Fourier transform of the Hamiltonian can be written as

$H({\bf k})=\begin{pmatrix}0&h({\bf k})\\\ h^{\dagger}({\bf
k})&0\end{pmatrix},h({\bf k})=\frac{t_{\sigma}}{2}\sum_{i=1}^{3}e^{ia{\bf
k}\cdot\hat{e}_{i}}\left(I+\sigma_{x}\cos 2\theta_{i}+\sigma_{y}\sin
2\theta_{i}\right),$ (S24)

where $\sigma_{x,y}$ are Pauli matrices in the chiral orbital basis
$(p_{-}/\sqrt{2},p_{+}/\sqrt{2})$, $p_{\mp}=p_{x}\mp ip_{y}$, $I$ is the
identity matrix, and we have defined $\theta_{1}=\pi/6$, $\theta_{2}=5\pi/6$,
and $\theta_{3}=3\pi/2$. Importantly, this Hamiltonian is symmetric under
time-reversal, inversion, and six-fold rotations.

The Hamiltonian in Eq. (S24) exhibits a quadratic band touching at ${\bf
k}=0$. Expanding to $\mathcal{O}(k^{2})$ we find

$h({\bf k})=\frac{3t_{\sigma}}{2}\begin{pmatrix}1&0\\\
0&1\end{pmatrix}+\frac{3t_{\sigma}a}{4}\begin{pmatrix}0&k_{+}\\\
-k_{-}&0\end{pmatrix}-\frac{3t_{\sigma}a^{2}}{16}\begin{pmatrix}2k^{2}&k_{-}^{2}\\\
k_{+}^{2}&2k^{2}\end{pmatrix}+\mathcal{O}(k^{3}).$ (S25)

After performing a unitary rotation in the sublattice space (two different
sublattices originate from non-equivalent sites $A$ and $B$ in the primitive
cell of the honeycomb lattice), we obtain

$\displaystyle H({\bf k})\to U^{\dagger}H({\bf
k})U\approx\frac{3t_{\sigma}}{2}\begin{pmatrix}1-\frac{(ak)^{2}}{4}&-\frac{(ak_{-})^{2}}{8}&0&\frac{ak_{+}}{2}\\\
-\frac{(ak_{+})^{2}}{8}&1-\frac{(ak)^{2}}{4}&-\frac{ak_{-}}{2}&0\\\
0&-\frac{ak_{+}}{2}&-1+\frac{(ak)^{2}}{4}&\frac{(ak_{-})^{2}}{8}\\\
\frac{ak_{-}}{2}&0&\frac{(ak_{+})^{2}}{8}&-1+\frac{(ak)^{2}}{4}\end{pmatrix},$

where $U=(I-i\tau_{y})/\sqrt{2}$ and $\tau_{y}$ is the Pauli matrix in the
sublattice space. Rewritten in this form, the Hamiltonian is block-diagonal:
the lower right (upper left) block describes quadratic band touching at energy
$-3t_{\sigma}/2$ ($3t_{\sigma}/2$), while the off-diagonal blocks describe
virtual hopping between the low- and high-energy bands. Accounting for such
hopping within second-order perturbation theory (?, ?), we end up with the
effective low-energy quadratic Hamiltonian describing just two lower bands
(which correspond to $n=3$ and $n=4$) that exhibit a QBTP:

$H_{\text{eff}}({\bf
k})\approx\frac{1}{4m}\left(\begin{array}[]{cc}k^{2}&k_{-}^{2}\\\
k_{+}^{2}&k^{2}\end{array}\right),$ (S26)

where the effective mass is given by $m^{-1}=3t_{\sigma}a^{2}/4$ and the
constant term $-(3t_{\sigma}/2)I$ has been dropped. The eigenenergies of this
Hamiltonian are given by ${\varepsilon}_{3}=0$ and
${\varepsilon}_{4}=k^{2}/2m$, i.e., we have a band touching between a
quadratically-dispersing band and a flat band. The eigenvectors equal

$|u_{\bf
k}^{3}\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}-k_{-}^{2}/k^{2}\\\
1\end{array}\right),\qquad|u_{\bf
k}^{4}\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}k_{-}^{2}/k^{2}\\\
1\end{array}\right),$ (S27)

resulting in the Berry connection of the form

${\bf A}_{{\bf k}}^{33}={\bf A}_{{\bf k}}^{44}=-{\bf A}_{{\bf k}}^{34}=-{\bf
A}_{{\bf
k}}^{43}=\frac{1}{k_{x}^{2}+k_{y}^{2}}\left(\begin{array}[]{c}-k_{y}\\\
k_{x}\end{array}\right).$ (S28)

The Berry phase around the QBTP is twice that around the LBTP, i.e., equals
$\pm 2\pi$. Since the Berry phase is defined only mod $2\pi$, such a
characteristic is not particularly useful in distinguishing a topologically-
trivial QBTP from a non-trivial one. Instead, as we discuss below, one can use
a gauge-independent integer winding number, provided certain symmetries are
present.

Repeating the same derivation as for the LBTP, we find that the population of
the $n=4$ band after scattering off the singularity (assuming the particles
are originally in the $n=3$ band) and the quantum distance are given by

$d^{2}_{\text{quad}}({\bf k},{\bf k}^{\prime})=|\langle u_{{\bf
k}^{\prime}}^{3}|u_{{\bf k}}^{4}\rangle|^{2}=\sin^{2}\theta.$ (S29)

It is worth discussing now certain important differences between the
geometries of the linear and quadratic band crossings. For a comprehensive
analysis of this and related questions we refer to Ref. (?), while we only
focus on the most important points.

In general, linear band touchings are always characterized by a Berry phase of
$\pi$ (mod $2\pi$) around them, and the maximal quantum distance between the
states encircling the point always equals $d^{\max}_{\text{Dirac}}=1$, meaning
that there are always orthogonal states no matter how close we are to the
singularity. The linear singularities are always topological in this sense and
no symmetries play any role in this conclusion.

Quadratic band touching points are quite different in many ways. Despite the
common belief, the Berry phase around the QBTP, in general, does not have to
be quantized to integers of $2\pi$. Related to this fact, the maximal quantum
distance between the states near the band touching point, in general, can take
an arbitrary value between 0 and 1, and so can the amplitude of the transition
probability oscillations $|\langle u_{{\bf k}^{\prime}}^{3}|u_{{\bf
k}}^{4}\rangle|^{2}$.

However, this arbitrariness is removed when certain symmetries are present in
the system. In particular, when time-reversal and rotational $C_{6}$
symmetries are present, which is the case in our system, the Berry phase
becomes quantized in units of $2\pi$, and the quantum distance can only take
values 0 or 1. We emphasize that the quantum distance is a more useful
quantity for characterizing the non-trivial topology of the band touching.
Indeed, while the Berry phase is defined only mod $2\pi$, its values of $\pm
2\pi$ can be easily changed to 0 by a simple gauge transformation. It implies
that the Berry phase cannot reliably distinguish between a trivial QBTP with
$d_{\max}=0$ and a non-trivial one with $d_{\max}=1$.

Equivalently, when time-reversal and rotational $C_{6}$ symmetries are
present, the non-trivial topology of the QBTP can be characterized by a well-
defined gauge-independent winding number (?). In general, the effective two-
band Hamiltonian can always be rewritten as a pseudo-spin in a pseudo-magnetic
field ${\bf B}({\bf k})$, $H_{\text{eff}}({\bf k})=B_{0}({\bf k})\cdot I+{\bf
B}({\bf k})\cdot{\boldsymbol{\sigma}}$, where ${\boldsymbol{\sigma}}$ is the
vector of Pauli matrices. If the aforementioned symmetries are present, one
can choose a basis such that $\left\\{{\bf B}({\bf
k})\cdot{\boldsymbol{\sigma}},\sigma_{z}\right\\}=0$, i.e., $B_{z}({\bf
k})=0$. In this case, the pseudo-magnetic field lies entirely in the
transverse pseudo-spin plane, where it traces out a great circle on the unit
Bloch sphere as ${\bf k}$ goes around the band touching point. In this sense,
the winding of ${\bf B}({\bf k})$ around the origin (singularity) is a
topological property since it cannot be changed smoothly. The winding number
counts how many times ${\bf B}({\bf k})$ encircles the origin and is well
defined as long as ${\bf B}({\bf k})\neq 0$ on a trajectory that encompasses
the singularity. It can be explicitly defined as (?, ?, ?)

$w=\frac{i}{2\pi}\oint d{\bf k}\nabla_{\bf k}\ln B_{-}({\bf k}),$ (S30)

where the contour of integration encircles the singularity and
$B_{-}=B_{x}-iB_{y}$. Near the QBTP in our system, $B_{0}({\bf k})=k^{2}/4m$,
${\bf B}({\bf k})=(1/4m)(k_{x}^{2}-k_{y}^{2},2k_{x}k_{y},0)^{T}$, and
$B_{-}({\bf k})=k_{-}^{2}/4m$, and so the winding number for the QBTP
$w_{\text{quad}}=2$. For comparison, the winding number for Dirac Hamiltonian
(S15) equals $w_{\text{Dirac}}=-1$.

The amplitude of the oscillations measured experimentally and shown in Fig. 4C
of the main text differs from 1 due to nonadiabaticity in the lattice
acceleration steps and, particularly, in the band mapping procedure.
Nevertheless, according to the discussion above, the non-zero amplitude and
two cycles of oscillations combined with our lattice symmetries unambiguously
indicate that the topological winding number around the QBTP we study is well-
defined and equals 2.