Determination of Chern numbers with a phase retrieval algorithm

TL;DR

This paper introduces a novel phase retrieval algorithm to measure Chern numbers in 2D topological insulators using ultracold bosonic atoms, enabling efficient and robust topological invariant detection via time-of-flight imaging.

Contribution

It presents a new experimental scheme combining phase retrieval with ultracold atoms to determine Chern numbers in topological insulators, improving measurement robustness.

Findings

Effective extraction of Chern numbers from time-of-flight images.

Demonstrated method's robustness with Harper-Hofstadter and Haldane models.

Applicable to multiband topological systems with bosonic atoms.

Abstract

Ultracold atoms in optical lattices form a clean quantum simulator platform which can be utilized to examine topological phenomena and test exotic topological materials. Here we propose an experimental scheme to measure the Chern numbers of two-dimensional multiband topological insulators with bosonic atoms. We show how to extract the topological invariants out of a sequence of time-of-flight images by applying a phase retrieval algorithm to matter waves. We illustrate advantages of using bosonic atoms as well as efficiency and robustness of the method with two prominent examples: the Harper-Hofstadter model with an arbitrary commensurate magnetic flux and the Haldane model on a brick-wall lattice.

Equations

$${\psi }_{\mathbf{k}}^{[n]}\left(\mathbf{r}\right)\propto \sum _{\ell ,\alpha }{e}^{i\mathbf{k}\cdot {\mathbf{r}}_{\ell \alpha }}{u}_{\alpha }^{[n]}\left(\mathbf{k}\right)w\left(\mathbf{r}-{\mathbf{r}}_{\ell \alpha }\right),$$

$$\mathcal{H}\left(\mathbf{k}\right){\mathbf{u}}^{[n]}\left(\mathbf{k}\right)=E^{[n]}\left(\mathbf{k}\right){\mathbf{u}}^{[n]}\left(\mathbf{k}\right),$$

$${\mathcal{A}}_{\mu }^{[n]}\left(\mathbf{k}\right)={\mathbf{u}}^{[n]†}\left(\mathbf{k}\right){\partial }_{\mu }{\mathbf{u}}^{[n]}\left(\mathbf{k}\right),$$

$$F_{xy}^{[n]}\left(\mathbf{k}\right)={\partial }_x{\mathcal{A}}_y^{[n]}\left(\mathbf{k}\right)-{\partial }_y{\mathcal{A}}_x^{[n]}\left(\mathbf{k}\right),$$

$$c_n=\frac{1}{2\pi i}{\int }_{\text{BZ}}{\mathrm{d}}^2\mathbf{k}{F}_{xy}^{[n]}\left(\mathbf{k}\right),$$

$${\sigma }_{xy}=-\left(e^2/h\right)\sum _n{c}_n,$$

$$I\left(\mathbf{r}\right)\propto |{\tilde{\psi }}_{\mathbf{k}}\left(\mathbf{q}\right){|}^2\propto {\left|\int \text{mbox}{d}^2\mathbf{r}{e}^{-i\mathbf{q}\cdot \mathbf{r}}{\psi }_{\mathbf{k}}\left(\mathbf{r}\right)\right|}^2,\mathbf{q}=\frac{m\mathbf{r}}{\hslash t_{\text{TOF}}},$$

$$\epsilon =\int {\mathrm{d}}^2\mathbf{q}{\left(\left|{\tilde{\psi }}^{\text{ ret}}\left(\mathbf{q}\right)\right|-\left|{\tilde{\psi }}_{\mathbf{k}}\left(\mathbf{q}\right)\right|\right)}^2,$$

$$u_{\alpha }\left(\mathbf{k}\right)=\frac{1}{n_{\text{cells}}}\sum _{\ell }{e}^{-i\mathbf{k}\cdot {\mathbf{r}}_{\ell \alpha }}\int {\mathrm{d}}^2\mathbf{r}{w}^{\ast }\left(\mathbf{r}-{\mathbf{r}}_{\ell \alpha }\right){\psi }_{\mathbf{k}}\left(\mathbf{r}\right).$$

$$|{\tilde{\psi }}_{\mathbf{k}}(\mathbf{q}){|}^2\propto |\tilde{w}\left(\mathbf{q}\right){|}^2{\left|\sum _{\ell ,\alpha }{e}^{i(\mathbf{k}-\mathbf{q})\cdot {\mathbf{r}}_{\ell \alpha }}{u}_{\alpha }\left(\mathbf{k}\right)\right|}^2,$$

$$\hat{\mathcal{H}}=-\sum _{m,n}\left(J{\hat{c}}_{m+1,n}^{†}{\hat{c}}_{m,n}+J^{\prime }{e}^{i2\pi \varphi m}{\hat{c}}_{m,n+1}^{†}{\hat{c}}_{m,n}+H.c.\right),$$

$$-J{e}^{i{k}_x}{u}_{\alpha +1}\left(\mathbf{k}\right)-2J^{\prime }$$

$$-J{e}^{-i{k}_x}{u}_{\alpha -1}\left(\mathbf{k}\right)=E\left(\mathbf{k}\right)u_{\alpha }\left(\mathbf{k}\right),$$

$$\hat{\mathcal{H}}=-J\sum _{⟨i,j⟩}{\hat{c}}_i^{†}{\hat{c}}_j-J^{\prime }\sum _{\text{llangle}i,j\text{rrangle}}{e}^{i{\theta }_{ij}}{\hat{c}}_i^{†}{\hat{c}}_j+\Delta \sum _i{ϵ}_i{\hat{c}}_i^{†}{\hat{c}}_i,$$

$${\mathcal{H}}_{11}=\Delta -2J^{\prime }(\cos \left(\theta +2k_x\right)+\cos \left(\theta -k_x-k_y\right)$$

$$+\cos \left(\theta -k_x+k_y\right)),$$

$${\mathcal{H}}_{12}={\mathcal{H}}_{21}^{\ast }=-J(2\cos k_x+e^{-i{k}_y}),$$

$${\mathcal{H}}_{22}=-\Delta -2J^{\prime }(\cos \left(\theta -2k_x\right)+\cos \left(\theta +k_x+k_y\right)$$

$$+\cos \left(\theta +k_x-k_y\right)).$$

$${\psi }_{\mathbf{k}}=\alpha {\psi }_{\mathbf{k}}^{[1]}+\beta {\psi }_{\mathbf{k}}^{[2]},{\left|\alpha \right|}^2+{\left|\beta \right|}^2=1,$$

$$\mathcal{H}\left(\mathbf{k}\right){\mathbf{u}}^{[n]}\left(\mathbf{k}\right)=E^{[n]}\left(\mathbf{k}\right){\mathbf{u}}^{[n]}\left(\mathbf{k}\right),$$

$${\mathbf{k}}_l$$

$$k_{l_{\mu }}$$

$$U_{\mu }^{[n]}\left({\mathbf{k}}_l\right):={\mathbf{u}}^{[n]†}\left({\mathbf{k}}_l\right){\mathbf{u}}^{[n]}\left({\mathbf{k}}_l+\hat{\mu }\right)/{\mathcal{N}}_{\mu }^{[n]}\left({\mathbf{k}}_l\right),$$

$${\tilde{F}}_{xy}^{[n]}\left({\mathbf{k}}_l\right):=\ln \frac{U_x^{[n]}\left({\mathbf{k}}_l\right)U_y^{[n]}\left({\mathbf{k}}_l+\hat{x}\right)}{U_x^{[n]}\left({\mathbf{k}}_l+\hat{y}\right)U_y^{[n]}\left({\mathbf{k}}_l\right)},$$

$$-\pi <\frac{1}{i}{\tilde{F}}_{xy}^{[n]}\left({\mathbf{k}}_l\right)\le \pi .$$

$$c_n=\frac{1}{2\pi i}\sum _l{\tilde{F}}_{xy}^{[n]}\left({\mathbf{k}}_l\right).$$

$$\epsilon =\int {\mathrm{d}}^2\mathbf{q}{\left(\left|{\tilde{\psi }}^{(i)}\left(\mathbf{q}\right)\right|-\left|{\tilde{\psi }}_{\mathbf{k}}\left(\mathbf{q}\right)\right|\right)}^2,$$

$${\mathbf{P}}_{w_1}{\psi }^{(i)}\left(\mathbf{r}\right):=\frac{1}{\mathcal{N}}\sum _{\ell ,\alpha }{v}_{\ell \alpha }^{(i)}w\left(\mathbf{r}-{\mathbf{r}}_{\ell \alpha }\right),$$

$$v_{\ell \alpha }^{(i)}=\int {\mathrm{d}}^2\mathbf{r}{w}^{\ast }\left(\mathbf{r}-{\mathbf{r}}_{\ell s}\right){\psi }^{(i)}\left(\mathbf{r}\right).$$

$${\mathbf{P}}_{w_2}\left({\mathbf{P}}_{w_1}{\psi }^{(i)}\left(\mathbf{r}\right)\right):=\sum _{\ell ,\alpha }\text{absolutevalue}{v_{\alpha }^{(i)}}_{\text{rms}}{e}^{i\mathrm{Arg}{v}_{\ell \alpha }^{(i)}}w\left(\mathbf{r}-{\mathbf{r}}_{\ell \alpha }\right),$$

$$\text{absolutevalue}{v_{\alpha }^{(i)}}_{\text{rms}}^2=\frac{1}{n_{\text{cells}}}\sum _{\ell }\text{absolutevalue}{v_{\ell \alpha }^{(i)}}^2$$

$${\mathbf{P}}_w={\mathbf{P}}_{w_2}\circ {\mathbf{P}}_{w_1},$$

Full text

I Introduction

Since Richard Feynman presented new perspectives of simulating physics Feynman (1982), there has been an outburst of works devoted to quantum simulators Buluta and Nori (2009); Hauke et al. (2012a); Cirac and Zoller (2012), which are relatively simple and controllable quantum systems that can experimentally emulate the behavior of other quantum systems or phenomena. A pronounced advantage of quantum simulators is most apparent when a targeted system is too difficult to handle for classical computers or when it is inaccessible experimentally.

Photonic devices Aspuru-Guzik and Walther (2012), trapped ions Blatt and Roos (2012) and ultracold atoms Bloch et al. (2008, 2012); Celi et al. (2017) are considered as the most promising quantum simulator platforms. In particular, ultracold atoms in optical lattices constitute clean feasible systems that are free from lattice defects, phonon vibrations and electron-electron interactions. As such, these systems seem to be especially well suited to mimic miscellaneous condensed matter phenomena Jaksch and Zoller (2005); Lewenstein et al. (2013); Dutta et al. (2015). By introducing fast periodic lattice modulations such as lattice shaking Eckardt et al. (2005) or laser-assisted tunneling Miyake et al. (2013) (for a review see Eckardt (2017)) it is possible to study classical magnetism Struck et al. (2011); Kosior and Sacha (2013) and create synthetic magnetic fields for neutral atoms Jaksch and Zoller (2003); Kolovsky (2011); Goldman et al. (2014); Goldman and Dalibard (2014); Celi et al. (2014) and successively design non-Abelian gauge potentials Tagliacozzo et al. (2013a); Kosior and Sacha (2014a), quantum simulators of lattice gauge theories Hauke et al. (2012b); Banerjee et al. (2012); Zohar et al. (2012); Banerjee et al. (2013); Tagliacozzo et al. (2013b); Zohar et al. (2013); Notarnicola et al. (2015); Kasper et al. (2016); Dutta et al. (2017); Zohar et al. (2017); González-Cuadra et al. (2017) and topologically non-trivial quantum systems Miyake et al. (2013); Jotzu et al. (2014); Aidelsburger et al. (2015); Kennedy et al. (2015); Stuhl et al. (2015); Mancini et al. (2015); Goldman et al. (2016); Kolovsky (2018).

The topologically protected edge conductivity in quantum Hall systems and in topological insulators is a consequence of topological properties of energy bands Halperin (1982); Hasan and Kane (2010); Qi and Zhang (2011). As in the celebrated Harper-Hofstadter model Harper (1955); Hofstadter (1976) and the Haldane model Haldane (1988) (for experiments in ultracold atoms see Miyake et al. (2013); Jotzu et al. (2014); Aidelsburger et al. (2015); Kennedy et al. (2015)), the energy bands are characterized by a non-zero value of topologically invariant Chern numbers. There are a few proposals how to measure the Chern numbers in a two dimensional (2D) ultracold quantum systems, including the center of mass motion Price and Cooper (2012); Dauphin and Goldman (2013); Aidelsburger et al. (2015); Dauphin et al. (2017) and direct time-of-flight (TOF) measurements with fermionic atoms Alba et al. (2011); Hauke et al. (2014); Fläschner et al. (2016) (see also other relevant works in strip geometries Wang et al. (2013); Schweizer et al. (2016); Lu et al. (2016); Mugel et al. (2017) and a very recent proposal on measuring Floquet topological invariants Ünal et al. (2018)).

In this paper, we propose an efficient method to determine Chern numbers of a 2D multi-band topological insulator in a series of standard TOF measurements with a single component Bose-Einstein Condensate (BEC) prepared in an optical lattice. We apply a phase retrieval algorithm Fienup (1978, 1982); Fienup and Wackerman (1986); Marchesini (2007); Kosior and Sacha (2014b) to matter waves in order to recover a small set of eigenstates that belong to the first Brillouin Zone (BZ). We illustrate robustness of the method with two important examples: the multiband Harper-Hofstadter model Harper (1955); Hofstadter (1976), with an arbitrary rational flux, and the Haldane model Haldane (1988) on the brick-wall lattice.

The paper is organized as follows. In Sec. II we present basic introduction to topological invariants of 2D Chern insulators and description of all elements of our method for determination of the Chern numbers. In Sec. III we show the main results of the numerical simulations demonstrating the application of the method. Section IV is devoted to an analysis of robustness of the method against experimental imperfections. We conclude in Sec. V.

II Method for determination of Chern numbers

We begin with a short introduction to Chern insulators and then we present all elements of the method for determination of Chern numbers in experiments with the help of a phase retrieval algorithm.

II.1 Topology of energy bands

Consider a general two-dimensional tight-binding model corresponding to a square optical lattice with the lattice spacings $a_{x}=a_{y}=1$. Assume that the Hamiltonian possesses discrete translational symmetries in the configuration space: $x\rightarrow x+q$ and $y\rightarrow y+1$ where $q$ is integer. In this case, a $q\times 1$ elementary cell has $q$ sublattice sites $\alpha=1,\dots,q$. Due to the translation symmetry, the system shows $q$ energy bands. An eigenstate belonging to the $n$-th band (where $n=1,\dots,q$) reads

[TABLE]

where $w\left(\bm{r}-\bm{r}_{\ell\alpha}\right)$ is the Wannier function localized at the site $\bm{r}_{\ell\alpha}=(\alpha,\ell)$ of the optical lattice, $\bm{k}=(k_{x},k_{y})$ is the system quasimomentum, where $k_{x}\in\left(-\pi/q,\pi/q\right]$ and $k_{y}\in\left(-\pi,\pi\right]$, and $u_{\alpha}^{[n]}=u_{\alpha+q}^{[n]}$ is a complex valued $q$-periodic function. Due to the translational symmetry of the system the full tight-binding Hamiltonian $\mathcal{H}$ can be written in a block diagonal form $\mathcal{H}=\bigoplus_{\bm{k}}\mathcal{H}\left(\bm{k}\right)$, where $\mathcal{H}\left(\bm{k}\right)$ are $q\times q$ blocks labeled by a quasimomentum $\bm{k}$ Bernevig and Hughes (2 13). The reduced Schrödinger equation

[TABLE]

where $\bm{u}^{[n]}\left(\bm{k}\right)=\left[u_{1}^{[n]}\cdots\;u_{q}^{[n]}\right]^{\top}$ is the normalized eigenvector, can be solved separately for each $\bm{k}$. The eigenenergies $E^{[n]}\left(\bm{k}\right)$ form a band. The geometry of energy bands can be described by the Berry connection $\mathcal{A}_{\mu}^{[n]}\left(\bm{k}\right)$ and Berry curvature $F_{xy}^{[n]}\left(\bm{k}\right)$ that read

[TABLE]

where $\mu=x,y$ denotes a direction in the quasimomentum space and $\partial_{\mu}=\partial/\partial k_{\mu}$ Bernevig and Hughes (2 13). Geometric features of energy bands can be related to topology - topological properties of the $n$-th band are characterized by the topologically invariant integer Chern number $c_{n}$, defined as an integral of the Berry curvature over the first BZ Bernevig and Hughes (2 13)

[TABLE]

The Chern numbers determine the Hall conductance of the system if fermions are loaded to the optical lattice. The total Hall conductance is the sum of conductances of all energy bands below the Fermi level and reads

[TABLE]

which is the famous Thouless-Kohmoto-Nightingale-den Nijs (TKKN) formula Thouless et al. (1982).

In practice, it is very efficient to calculate the Chern number using the Fukui-Hatsugai-Suzuki (FHS) method Fukui et al. (2005) which, rather then a crude discretization of (5), exploits the lattice gauge theory formalism by defining the Berry connection on the coarsely discretized BZ (see Appendix A for details).

In the following we show that applying the FHS approach and a phase retrieval algorithm Fienup (1978, 1982); Fienup and Wackerman (1986); Marchesini (2007); Kosior and Sacha (2014b) we can reconstruct Chern numbers from a series of time-of-flight experiments with a single component BEC.

II.2 Preparation of initial eigenstates

If we knew all eigenstates of a given energy band of the Hamiltonian, then Eq. (5) would allow us to obtain the Chern number characterizing the band. We will show that when a BEC in the optical lattice is prepared in a certain eigenstate, measurement of the density of atoms after TOF and application of a phase retrieval algorithm allow us to reconstruct the wavefunction completely. Performing the same experiments but with a BEC in different eigenstates of the band provides sufficient information to determine the Chern number of the band. (See Sec.IV A for the analysis of the BZ meshing size.) In this subsection we discuss the first element of the method, i.e. the preparation of a BEC in different eigenstates of an energy band cFi .

To prepare a BEC in an eigenstate corresponding to a topologically non-trivial energy band, one usually starts an experimental sequence with loading a BEC into the ground state of a 2D optical lattice with trivial topology Aidelsburger et al. (2015). The ground state can be well-approximated by a Bloch wave (1) with a quasimomentum $\bm{k}_{in}$ that minimizes the dispersion relation. By turning on artificial gauge fields, the system is then driven into a regime of non-trivial topology of energy bands which are characterized by non-zero values of the Chern numbers (5). However, while switching from trivial to non-trivial topology, a quantum phase transition takes place which is accompanied by closing a gap between a neighboring band at distinct quasimomenta $\bm{k}_{D}\in\mathcal{D}$ (the set of Dirac points) Bernevig and Hughes (2 13). If $\bm{k}_{in}\approx\bm{k}_{D}$, in order to avoid population of another band, before we change parameters of the system across the topological quantum phase transition, we have to apply a weak constant force $\bm{F}_{1}$ for a suitable period of time $\Delta t_{1}$ so that the system is transferred to some auxiliary quasi-momentum $\bm{k}_{aux}=\bm{k}_{in}+\frac{\Delta t_{1}}{\hbar}\bm{F}_{1}\neq\bm{k}_{D}$ (see Fig. 1). Then, slow change of parameters of the system across the topological phase transition does not lead to population of another band if it is done on a time scale longer than the scale given by the inverse of the energy gap corresponding to $\bm{k}_{aux}$. Once we are in the topological phase, we can apply another weak force $\bm{F}_{2}$ which allows us to transfer the system to any quasi-momentum $\bm{k}=\bm{k}_{aux}+\frac{\Delta t_{2}}{\hbar}\bm{F}_{2}$ we need. In Sec. II.3 we show how to recover full information about an eigenstate of the Bose system corresponding to a given quasi-momentum $\bm{k}$ in the measurement of the atomic density after TOF. Following this experimental sequence, we can scan the whole first BZ in separate experimental realizations and obtain sufficient information about the system which allows one to determine the Chern numbers by means of the FHS approach. In the presented experimental scheme we argue that using bosonic atoms it is possible to switch to the non-trivial topology almost adiabatically by avoiding band touching points. Nevertheless, in Sec. IV.2 we present numerical studies of the influence of excitations to other bands on the determination of the Chern numbers.

II.3 Phase retrieval after TOF

In this section we review and adapt a method Kosior and Sacha (2014b) which allows one to reconstruct a BEC wavefunction out of a standard time-of-flight image after being processed with a phase retrieval algorithm Fienup (1978, 1982); Fienup and Wackerman (1986); Marchesini (2007).

A time-of-flight image shows the spatial density distribution $I\left(\bm{r}\right)$ of atoms after a time period $t_{TOF}$ of a free expansion that follows a sudden turning off an optical lattice and external trapping potentials. In the far field limit, $I\left(\bm{r}\right)$ is proportional to the initial distribution of atoms in the momentum space if we may neglect interaction between particles during the expansion of the atomic cloud Pedri et al. (2001); Gerbier et al. (2008); Deuar (2016); cSe

[TABLE]

where $\bm{k}$ is the initial quasimomentum, $\psi_{\bm{k}}\left(\bm{r}\right)$ and $\tilde{\psi}_{\bm{k}}\left(\bm{q}\right)$ are the representations of the initial condensate wave function in the real and reciprocal spaces, and $m$ is the atomic mass. A measurement of the atomic density reveals $|\tilde{\psi}_{\bm{k}}\left(\bm{q}\right)|^{2}$ at discrete points in the $\bm{q}$ space. If we knew not only the density but also the phase of $\tilde{\psi}_{\bm{k}}\left(\bm{q}\right)$ we would be able to obtain the wavefunction $\psi_{\bm{k}}\left(\bm{r}\right)$ by means of the inverse discrete Fourier transform. However, even without the knowledge of the phase, the task is not hopeless if we have some additional information about the system. Ultra-cold atoms are always prepared in a trap, i.e. the system always occupies finite area in the configuration space. If the support $S$ of $\psi_{\bm{k}}\left(\bm{r}\right)$ (area where $\psi_{\bm{k}}\left(\bm{r}\right)\neq 0$) and the modulus $|\tilde{\psi}_{\bm{k}}\left(\bm{q}\right)|$ are known, one can employ an iterative phase retrieval algorithm Fienup (1978, 1982); Fienup and Wackerman (1986); Marchesini (2007); Kosior and Sacha (2014b). Let us stress here that the presence of an external trap is indispensable but its shape is not important as long as the trap size is significantly larger than the lattice spacing so that the quasimomentum is a good quantum number. In the present article we consider ultra-cold atoms in optical lattices and in the presence of an external hard wall potential but the phase retrieval algorithm can be applied to other trapping potentials and lattice geometries. For example, in Ref. Kosior and Sacha (2014b) a 2D triangular lattice and a harmonic trapping potential are analyzed within a Thomas-Fermi approximation, where the external potential modifies the envelope of the wavefunction only.

The phase retrieval algorithm seeks for the intersection of two sets of functions: a set of functions with a given support $S$ in the position space and a set of functions with a given modulus $|\tilde{\psi}_{\bm{k}}\left(\bm{q}\right)|$ in the reciprocal space. Let $\psi^{(i)}\left(\bm{r}\right)$ be an approximation of the desired solution at $i$-th iteration of the phase retrieval algorithm. The algorithm starts with a random, complex-valued $\psi^{(0)}\left(\bm{r}\right)$ that satisfies the support constraint $\psi^{(0)}\left(\bm{r}\right)=0$ for $\bm{r}~{}\notin~{}S$. In the simplest version of the algorithm Fienup (1978), the following operations are performed at each iteration:

1. (i)

The Fourier transform is performed on $\psi^{(i)}\left(\bm{r}\right)$, resulting in $|\tilde{\psi}^{(i)}\left(\bm{q}\right)|e^{i\phi^{(i)}\left(\bm{q}\right)}$. 2. (ii)

$|\tilde{\psi}^{(i)}\left(\bm{q}\right)|$ is substituted with the true $|\tilde{\psi}_{\bm{k}}\left(\bm{q}\right)|$ which is obtained in an experiment after TOF. 3. (iii)

Inverse Fourier transform is applied which gives $\psi^{(i+1)}\left(\bm{r}\right)$, not necessarily satisfying the support constraint. 4. (iv)

The support constraint is imposed on $\psi^{(i+1)}\left(\bm{r}\right)$ by setting $\psi^{(i+1)}\left(\bm{r}\right)=0$ for every $\bm{r}\notin S$.

Convergence of the algorithm is tracked by the error measure defined as

[TABLE]

where $\tilde{\psi}^{\text{ret}}_{\bm{k}}\left(\bm{q}\right)$ is a retrieved function and $\left|\tilde{\psi}_{\bm{k}}\left(\bm{q}\right)\right|^{2}$ is the measured probability distribution. The presented simplest version of the algorithm guarantees a decrease of $\varepsilon$ in each iteration. Unfortunately, once it reaches a local minimum of $\varepsilon$, it cannot proceed further. There are modifications of the phase retrieval methods which allow for the much faster convergence to a desired solution Fienup (1982); Fienup and Wackerman (1986); Marchesini (2007). Moreover, to increase the rate of the convergence one can use any extra information about $\psi_{\bm{k}}\left(\bm{r}\right)$, e.g., a preliminary in-situ measurement of $\left|\psi_{\bm{k}}\left(\bm{r}\right)\right|$ or its theoretical estimation Kosior and Sacha (2014b). In our case, we speed up the convergence by exploiting information about geometry of an optical lattice only, i.e. we do not assume anything about the parameters of the Hamiltonian, see Appendix B for all details.

Once $\psi_{\bm{k}}\left(\bm{r}\right)$, Eq. (1), is successfully recovered, in order to extract the coefficient vector $\bm{u}\left(\bm{k}\right)$ one has to project $\psi_{\bm{k}}\left(\bm{r}\right)$ on the orthonormal basis of the Wannier functions. To minimize the numerical error one might additionally average each $u_{\alpha}$ component over lattice sites $\ell=1,\ldots,n_{\text{cells}}$:

[TABLE]

The Wannier functions $w(\bm{r})$ can be well approximated by Gaussian functions with the width $\sigma_{w}=\hbar t_{\text{TOF}}/\left(m\sigma_{\tilde{w}}\right)$ which can be obtained from the wide envelope of the measured density profile

[TABLE]

where $\tilde{w}\left(\bm{q}\right)$ is the Fourier transform of $w(\bm{r})$ and $\sigma_{\tilde{w}}$ is the width of $\tilde{w}\left(\bm{q}\right)$.

II.4 Calculation of the Chern number

In order to determine the Chern number we propose a series of experiments with a BEC in an optical lattice. In each experiment, one prepares a BEC in an eigenstate with a different quasimomentum $\bm{k}$ from the first BZ and retrieves a column complex-valued vector $\bm{u}\left(\bm{k}\right)$, Eq. (2), using the phase retrieval algorithms (see Sec. II.3). To obtain the Chern number we apply a highly effective FHS method Fukui et al. (2005) which allows us to calculate the Chern number with the help of a few eigenvectors $\psi_{\bm{k}}(\bm{r})$ only, i.e. the coarsely discretized BZ. It is possible due to the fact that the FHS algorithm is based on a gauge-invariant lattice gauge theory formulation. (See Appendix A for a quick revision of the FHS algorithm.) In Sec. III we demonstrate the method of the determination of the Chern numbers simulating experimental data for two examples: Harper-Hofstadter and Haldane models.

III Numerical simulations

The proposed experimental scheme of detecting Chern numbers applies to a general tight-binding Hamiltonian in a two-dimensional space. In this section we illustrate application of the scheme with two examples: the Harper-Hofstadter model with an arbitrary rational flux Harper (1955); Hofstadter (1976) and the Haldane model Haldane (1988) on a brick-wall lattice (for experiments in ultracold atoms see Miyake et al. (2013); Jotzu et al. (2014); Aidelsburger et al. (2015); Kennedy et al. (2015)). In the case of the Harper-Hofstadter model we show that a large number of bands is not the limitation of our method. With the help of the Haldane model we demonstrate that our scheme allows one to reconstruct the phase diagram of the system.

III.1 The Harper-Hofstadter model

Consider bosonic atoms in a square two-dimensional lattice, in $XY$-plane, with a unit lattice spacing subjected to uniform artificial magnetic field $\bm{B}=\left(0,0,B\right)$. The nearest-neighbor-hopping Hamiltonian of an atom in the Landau gauge with the vector potential $\bm{A}=(0,Bx,0)$ takes the following form

[TABLE]

where $\hat{c}^{\dagger}_{m,n},\hat{c}_{m,n}$ are the bosonic particle creation and annihilation operators corresponding to a lattice site $(m,n)$. $J,\;J^{\prime}$ are tunneling amplitudes and $\phi=B/h$ is a dimensionless flux. Due to the presence of the magnetic field, particles tunneling along $y$ acquire the Peierls phase factor $e^{i2\pi\phi m}$ Peierls (1933). The presence of the magnetic field, in principle, breaks discrete space-translation symmetry of the lattice. However, if the flux is a rational number, $\phi=p/q$ where $p$ and $q$ are coprime integers, the translational symmetry is restored but with the spatial period $q$ times longer than the lattice constant. Therefore, an effective magnetic $q\times 1$ elementary cell consists of $q$ lattice sites, and the first BZ is the rectangle $2\pi/q\times 2\pi$ in the quasi-momentum space. After rewriting the Hamiltonian (11) in the Fourier space, the reduced Schrödinger equation (2) takes the following form:

[TABLE]

where $\alpha=1,\ldots,q$.

Let us focus on the reconstruction of the lowest band Chern number for $q=3$ and $q=5$ band models, as depicted in Fig. 2. We choose a finite optical lattice consisting of $7\times 7$ effective magnetic elementary cells (which corresponds to $7\times 21$ or $7\times 35$ lattice sites for $q=3$ and $q=5$ bands respectively). In principle, the measurement of $|\tilde{\psi}_{\bm{k}}(\bm{q})|$ and performing the phase retrieval algorithm allows us to obtain the full information about the eigenstate $\psi_{\bm{k}}(\bm{r})$ and successively recover the Chern number (see Sec. II.3 - II.4). However, the phase retrieval algorithm is known to occasionally get stuck at local minima. Therefore, for every $\bm{k}$ we repeat the algorithm, each time starting from different randomly generated initial state.

All retrieved eigenstates can be sorted by their error $\varepsilon$, Eq. (8), as shown in Fig. 3(upper panel). It is evident that about the 90% of the best phase retrieval runs converge to functions with approximately the same error $\varepsilon\approx 10^{-6}-10^{-7}$, while the errors of the last 5-10% trails are larger by a few orders of magnitude.

For each quasimomentum $\bm{k}$ we select a random representative out of 90 phase retrieval algorithm runs and calculate the Chern number with the FHS method. We repeat the process $10^{3}$ times and successively average the data. (Let us stress that this repetition is a data processing post measurement only.) As we illustrate in Fig. 3(lower panel), after rejecting the worst phase retrieval trails we are always able to recover the Chern numbers $c_{1}=-1$ with a perfect accuracy. Note that without any rejections, for a 8$\times$8 BZ mesh we obtain $c_{1}=-0.97(17)$ for $\phi=1/5$ and $c_{1}=-0.98(15)$ for $\phi=1/3$. Moreover, in Fig. 3 we show that a much harsher discretization of the first BZ (4$\times$4 mesh is already sufficient to correctly recover the Chern number.

III.2 Haldane model on a brick-wall lattice

The brick-wall structure consists of two interpenetrating square lattices $A$ and $B$, see Fig. 4. We assume real tunneling amplitudes $J$ between nearest neighboring lattice sites and complex tunneling amplitudes $J^{\prime}e^{\pm i\theta}$ between next-nearest neighboring sites. The model is topologically equivalent to the Haldane model on a honeycomb lattice Haldane (1988). The Hamiltonian of the system reads

[TABLE]

where $i,j$ are indices of the lattice sites, $\left\langle i,j\right\rangle$ denotes pairs of nearest neighbors, \left\llangle i,j\right\rrangle pairs of next-nearest neighbors, $\theta_{ij}=\pm\theta$ where the sign depends on the direction of the tunneling, $\Delta$ introduces the energy offset between the $A$ and $B$ sublattices because $\epsilon_{i}=1$ for $i\in A$, $\epsilon_{i}=-1$ for $i\in B$ (see Fig. 4). Complex values of the tunneling amplitudes break the time-reversal symmetry while the energy offset breaks the parity symmetry. Switching to the reciprocal space we can write the Hamiltonian in a block diagonal form. Each block is a $2\times 2$ matrix $\mathcal{H}(\bm{k})$ whose elements take the form

[TABLE]

An identical procedure as in the case of the Harper-Hofstadter model leads to a successful retrieval of the Chern number of the lowest band. This allows us to obtain the topological phase diagram of the Haldane model, see Fig. 5. The discretization of the first BZ corresponds to the $6\times 6$ mesh. For each of the eigenstates we assume that we know the support of $\psi_{\bm{k}}(\bm{r})$ and the modulus $|\tilde{\psi}_{\bm{k}}(\bm{q})|$ and perform the phase retrieval procedure 90 times. Each application of the algorithm starts with randomly chosen phases of an eigenstate and consists of 350 iterations. We may now select a number of the best results, based on their error $\varepsilon$, Eq. (8), and make statistics on the retrieved Chern numbers, as in Sec. III.1. Taking all results, including those that did not converge to a global solution, we obtain a topological phase diagram in Fig. 5 (upper panel) which only qualitatively represents a structure predicted by Haldane Haldane (1988). However, selecting 50% of the best results yields a perfect recovery of the Haldane model phase diagram, shown in Fig.5 (lower panel).

IV Robustness

In this section we investigate the influence of possible experimental imperfections on values of the retrieved Chern numbers. As an example we choose the Harper-Hofstadter Hamiltonian (11) with the flux $\phi=1/3$ and the finite lattice consisting of $7\times 7$ elementary magnetic cells ($7\times 21$ lattice sites). All presented quantities are averaged over 90 phase retrieval runs which correspond to different randomly chosen initial states. Percentage of discarded worst (according to error $\varepsilon$, Eq. (8)) retrieval results is either 10% or 90%. The error bars are the standard deviations of the averaged values.

IV.1 Number of points chosen in the first Brillouin Zone

We have tested how densely one has to probe the first BZ in order to get the proper value of the Chern number $c_{1}$ corresponding to the lowest energy band in Fig. 2 (left panel). Figure 6(a) indicates that it is sufficient to perform the $4\times 4$ mesh discretization of the Brillouin zone and the retrieved Chern number is correct. It also demonstrates how powerful the FHS method is. In order to make sure that a Chern number is retrieved correctly, an experiment should be performed again with different discretization of the BZ.

IV.2 Excitations to the second band

Experimental preparation of an eigenstate from the lowest energy band is usually not perfect and contributions from the higher bands can be expected. In this subsection we analyze contamination of eigenstates of the lowest (first) band $\psi_{\bm{k}}^{[1]}$ by eigenstates from the second band $\psi_{\bm{k}}^{[2]}$,

[TABLE]

and its influence on the determination of the Chern number $c_{1}$, which would estimate the worst case scenario for the Landau-Zener transition, see Sec. II.2.

We have applied our method for different populations $|\beta|^{2}$ of the second band and the results are presented in Fig. 6(b). We conclude that $\left|\beta\right|^{2}\lesssim 0.12$ allows for the correct retrieval of the Chern number $c_{1}$ in case of $8\times 8$ mesh. A similar, symmetric result applies to the Chern number of the second band: to obtain successfully $c_{2}$ we require $\left|\beta\right|^{2}\gtrsim 0.86$. If one takes only $4\times 4$ mesh, $\beta$ must satisfy $\left|\beta\right|^{2}\lesssim 0.02$ to recover the Chern number of the lowest band. Note that the mesh size in FHS method must be increased with the absolute value of the Chern number Fukui et al. (2005), and therefore in the case of the $4\times 4$ mesh it is not enough to recover a correct Chern number $c_{2}=2$. Although we find that the higher mesh gives a better critical $\left|\beta\right|^{2}$, at some point the undesired occupancies of other bands will always spoil the results. Therefore, in Sec. II.2 we propose a method to minimize the excitations to higher bands.

IV.3 Background noise

In the experiment, background noise will affect the atomic density measurements. Let us define the signal strength $A$ as the average value of $\left|\tilde{\psi}_{\bm{k}}(\bm{q})\right|^{2}$ calculated in the first BZ. The signal-to-noise ratio reads $\text{SNR}=A/\sigma_{n}$, where $\sigma_{n}$ is the standard deviation of Gaussian white noise whose absolute values are added to each point $\bm{q}$ of the atomic density image. The results of the retrieved Chern number versus SNR are presented in Fig. 6(c). The minimal SNR that allows for the successful retrieval of the Chern number is about 5.5 for a $8\times 8$ mesh after discarding about 90 percent of the worst retrievals. It is also important to note that experimental noise can be reduced either by repeating the experiment and averaging the recorded density profiles over separate realizations, or by applying noise removal algorithms V.R. et al. (2010); Lenzen et al. (2013); Niu et al. (2018).

IV.4 Resolution of experimental imaging system

In order to check how the results are sensitive to the resolution of the imaging system, we convolve the original atomic density after TOF, $|\psi_{\bm{k}}(\bm{q})|^{2}$, with the Gaussian profile of width $\sigma_{r}$. In Fig. 6(d) we can see how the average value of the Chern number $c_{1}$ depends on the ratio $\sigma_{r}/\sigma_{\text{peak}}$, where $\sigma_{\text{peak}}$ is the width of the Gaussian fit to the highest Bragg peak that can be observed in the atomic density, $|\psi_{\bm{k}}(\bm{q})|^{2}$, after TOF. The minimal resolution that guarantees the correct value of the Chern number is $\sigma_{r}/\sigma_{\text{peak}}\approx 0.4$, which also requires a $8\times 8$ mesh after discarding about 90% of the worst retrieval results.

V Summary

We have proposed a method for determination of the topological invariants of two-dimensional Chern insulators with the help of ultra-cold bosonic atoms in optical lattice potentials. The method relies on a sequence of experiments where a Bose-Einstein condensate is prepared in different eigenstates of a given energy band. In each experiment, an atomic density after time-of-flight is measured. Because the time-of-flight is actually the Fourier transform of the initial condensate wavefunction of atoms prepared in a finite optical lattice, a phase retrieval algorithm can be applied in order to obtain the phase of the wavefunction. The full knowledge of eigenstates of a given band allows one to calculate the Chern number characterizing the band.

We illustrate the application of the method with two examples: the Harper-Hofstadter model and the Haldane model on a brick wall lattice. It turns out that it is sufficient to retrieve a small number of eigenstates of a given band, i.e. to discretize coarsely the first Brillouin zone, in order to determine the Chern number. An experimental sequence that allows one to avoid population of neighboring bands, during the preparation of the system in a topological phase, is presented. We also analyze robustness of the method and its resistance to experimental imperfections.

ACKNOWLEDGMENTS

Support of Talenty and Grand Scholarship Programs (T.S.), the National Science Centre, Poland via Projects No. 2016/21/B/ST2/01095 (T.S.), No. 2016/21/B/ST2/01086 (A.K.) and QuantERA programme No. 2017/25/Z/ST2/03027 (K.S.) is acknowledged.

Appendix A Fukui-Hatsugai-Suzuki Method

Assume the 2D system on a square lattice that is invariant under discrete space translations $x\rightarrow x+q_{x}$ and $y\rightarrow y+q_{y}$, where $q_{x}$, $q_{y}$ are integer multiples of the lattice constant $a=1$. Hence, the system can be described completely by a $q_{x}q_{y}\times q_{x}q_{y}$ Hamiltonian matrix $\mathcal{H}\left(\bm{k}\right)$ in a reduced Brillouin zone $\bm{k}\in\left(-\pi/q_{x},\pi/q_{x}\right]\times\left(-\pi/q_{x},\pi/q_{x}\right]$. Assume that, for each $\bm{k}$, the Hamiltonian $\mathcal{H}\left(\bm{k}\right)$ has non-degenerate eigenvalues. Then, the solutions of the Schrödinger equation

[TABLE]

describe separate energy bands labeled by $n=1,\ldots,q_{x}q_{y}$. Let us take a set of discrete points $\bm{k}_{l}$ $\left(l=1,\dots,N_{x}N_{y}\right)$ in the first BZ

[TABLE]

with

[TABLE]

where $\mu=x,y$. We will call $\hat{\mu}$ the vector of the length $\delta k_{\mu}=2\pi/\left(q_{\mu}N_{\mu}\right)$ in the direction $\mu$. The $U(1)$ linking variables of the $n$-th band are defined as

[TABLE]

with $\mathcal{N}^{[n]}_{\mu}\left(\bm{k}_{l}\right)=\left|\bm{u}^{[n]\dagger}\left(\bm{k}_{l}\right)\bm{u}^{[n]}\left(\bm{k}_{l}+\hat{\mu}\right)\right|.$

The field strength $\tilde{F}^{[n]}_{xy}\left(\bm{k}_{l}\right)$ takes a manifestly gauge-invariant form

[TABLE]

Finally, the Chern number reads

[TABLE]

Even for coarsely discretized BZ’s this algorithm gives accurate values of the Chern numbers (see Sec. IV.1 or ref. Fukui et al. (2005)).

Appendix B Phase retrieval algorithm and its optimization

Phase retrieval algorithms iteratively seek for a solution $\psi\left(\bm{r}\right)$ in the object space, provided the modulus of its Fourier transform $\left|\tilde{\psi}_{\bm{k}}\left(\bm{k}\right)\right|$ and support S (area where $\psi\left(\bm{r}\right)\neq 0$) are known. The simplest version of the algorithm, called error reduction (ER), is described in Sec. II.3. Fienup proves Fienup (1978) that at each iteration, the retrieval error, defined as

[TABLE]

decreases. Stagnation of this algorithm in local minima is, however, likely to occur Fienup (1982); Marchesini (2007). Several approaches have been proposed to solve this problem Marchesini (2007). One example is the hybrid input-output (HIO) algorithm based on nonlinear feedback control theory Fienup (1982). It is very similar to the ER algorithm, the only change is the step (iv) described in Sec. II.3. The part of $\psi^{(i)}\left(\bm{r}\right)$ that lies outside the support is not set to zero but instead to $\left(1-\eta\bm{P}_{m}\right)\psi^{(i)}\left(\bm{r}\right)$, where the operator $\bm{P}_{m}$ (described in steps (i)-(ii) in Sec. II.3) is the projection on the set of functions with the modulus $\left|\tilde{\psi}_{\bm{k}}\left(\bm{k}\right)\right|$ and $0<\eta<1$ is the feedback parameter, usually set to $0.7-0.9$. In most cases, a combination of the HIO and ER methods, e.g. 20 iterations of HIO and 1 iteration of ER algorithms, repeated in cycle, gives the best results. Since the HIO method does not guarantee the decrease of the error $\varepsilon$, the last few (30-50) iterations, should consist of the pure ER algorithm.

Support

If we want to recover $N$ complex numbers $\psi_{\bm{k}}\left(\bm{r}\right)$ within support, we need at least $2N$ real numbers $\left|\tilde{\psi}_{\bm{k}}\left(\bm{q}\right)\right|$. This gives a constraint on the area of the support which must not be less than $50\%$ the area of the whole table of $\psi_{\bm{k}}\left(\bm{r}\right)$. In our case, the support occupies only $22.5\%$ of the whole table which increases the rate of convergence. If the support is symmetric with respect to rotation by 180 degrees around some point $\bm{r}_{0}$ in space (e.g., the support is a rectangle or a circle), the fact that $\psi_{\bm{k}}\left(\bm{r}-\bm{r}_{0}\right)$ and $\psi^{*}_{\bm{k}}\left(-\left(\bm{r}-\bm{r}_{0}\right)\right)$ have the same modulus of the Fourier transform causes an ambiguity. The algorithm will converge to any of the two solutions with equal probability and in some cases it will stagnate at their superposition Fienup and Wackerman (1986). The only other nonuniqness can appear if and only if $\psi_{\bm{k}}(\bm{r})$ can be written as a convolution of two non-central symmetric functions Barakat and Newsam (1984). Therefore, in our simulations we choose a trapezoidal support with the ratio 4/5 of its bases which corresponds to a hard-wall box potential of this shape. Fluctuations of the size of an atomic cloud in a trap result in changes of the width of the Bragg peaks in the momentum distribution. The latter are not dangerous in the determination of Chern numbers. We also stress that if the size of the cloud is fluctuating one must set a support that is slightly larger than the average size of the cloud. This way one does not unintentionally ”cut” the solution in real space.

Optimization

If additional information about $\psi\left(\bm{r}\right)$ is known, it can be used to speed up the algorithm convergence. For example if the geometry of an optical lattice and the number of lattice sites can be estimated in the experiment, we know all information about an eigenstate of the system presented in Eq. (1) except the factors $e^{i\bm{k}\cdot\bm{r}_{\ell\alpha}}u_{\alpha}$. We use this information as follows.

Define the projection $\bm{P}_{w_{1}}\psi^{(i)}$ of a current estimate of the desired solution on the Wannier state basis,

[TABLE]

where $\ell=1,\ldots,n_{\text{cells}}$ is the index of an elementary cell, $\alpha~{}=~{}1,\ldots,q$ is the index of a lattice site within an elementary cell, $\mathcal{N}$ is the normalization factor and

[TABLE]

If $\psi^{(i)}\left(\bm{r}\right)$ is identical with the desired solution, then $v_{\ell\alpha}^{(i)}\equiv e^{i\bm{k}\cdot\bm{r}_{\ell\alpha}}u_{\alpha}$, hence $|v_{\ell\alpha}^{(i)}|$ should not depend on $\ell$ and we impose this condition in the iterative process. We define the next projection $\bm{P}_{w_{2}}$,

[TABLE]

where

[TABLE]

is the mean occupation of the sublattice site $\alpha$. This operation ensures that occupations of the same sublattice sites in all elementary cells are the same (see Fig. 7 for clarification). The complete projection

[TABLE]

is performed every 3 iterations of the phase retrieval algorithm. The effect of our optimization is clearly visible in Fig. 8 — the final error (19) is about 4 orders of magnitude smaller than without the optimization (see also comprehensive phase retrieval software libraries Chandra et al. (2017)).