The dynamics of weakly interacting bosons in optical lattices with flux

The dynamics of weakly interacting bosons in optical lattices with flux

Ana Hudomal Scientific Computing Laboratory, Center for the Study of Complex Systems, Institute of Physics Belgrade, University of Belgrade, Serbia    Ivana Vasić Scientific Computing Laboratory, Center for the Study of Complex Systems, Institute of Physics Belgrade, University of Belgrade, Serbia    Hrvoje Buljan Department of Physics, Faculty of Science, University of Zagreb, Croatia    Walter Hofstetter Institut für Theoretische Physik, Johann Wolfgang Goethe-Universität, Frankfurt am Main, Germany    Antun Balaž Scientific Computing Laboratory, Center for the Study of Complex Systems, Institute of Physics Belgrade, University of Belgrade, Serbia
Abstract

Realization of strong synthetic magnetic fields in driven optical lattices has enabled implementation of topological bands in cold-atom setups. A milestone has been reached by a recent measurement of a finite Chern number based on the dynamics of incoherent bosonic atoms. The measurements of the quantum Hall effect in semiconductors are related to the Chern-number measurement in a cold-atom setup, however, the design and complexity of the two types of measurements are quite different. Motivated by these recent developments, we investigate the dynamics of weakly interacting incoherent bosons in a two-dimensional driven optical lattice exposed to an external force, which provides a direct probe of the Chern number. We consider a realistic driving protocol in the regime of high driving frequency and focus on the role of weak repulsive interactions. We find that interactions lead to the redistribution of atoms over topological bands both through the conversion of interaction energy into kinetic energy during the expansion of the atomic cloud and due to an additional heating. Remarkably, we observe that the moderate atomic repulsion facilitates the measurement by flattening the distribution of atoms in the quasi-momentum space. Our results also show that weak interactions can suppress the contribution of some higher-order non-topological terms in favor of the topological part of the effective model.

I Introduction

Ultracold atoms in optical lattices provide a perfect platform for quantum simulations of various condensed-matter phenomena Bloch et al. (2008). Yet, since charge-neutral atoms do not feel the Lorentz force, a big challenge in this field was realization of synthetic magnetic fields. After years of efforts, artificial gauge potentials for neutral atoms were implemented by exploiting atomic coupling to a suitable configuration of external lasers Lin et al. (2009); Dalibard et al. (2011). These techniques were further extended to optical lattices, leading to the realization of strong, synthetic, magnetic fields. As a result, important condensed-matter models - the Harper-Hofstadter Hofstadter (1976) and the Haldane model Haldane (1988) - are nowadays available in cold-atom setups Aidelsburger et al. (2013); Miyake et al. (2013); Jotzu et al. (2014). The key property of these models is their non-trivial topological content. In the seminal TKNN paper Thouless et al. (1982) it was shown that the quantization of the Hall conductivity observed in the integer Hall effect can be directly related to the topological index of the microscopic model - the Chern number.

Cold-atom realizations of topological models exploit periodic driving, either through laser-assisted tunneling Aidelsburger et al. (2013); Miyake et al. (2013) or by lattice shaking Jotzu et al. (2014). Using Floquet theory Floquet (1883); Grifoni and Hänggi (1998), a periodically driven system can be related to the time-independent effective Hamiltonian that corresponds to a relevant condensed-matter system. The mapping is known as Floquet engineering and its important features in the context of optical lattices are discussed in Refs. Goldman and Dalibard (2014); Goldman et al. (2015); Eckardt and Anisimovas (2015); Aidelsburger et al. (2017); Eckardt (2017); Cooper et al. (2018); Sun and Eckardt (2018); Fujiwara et al. (2018). Because of important differences of cold-atom setups and their condensed-matter counterparts, new quench protocols for probing topological features were proposed Price and Cooper (2012); Dauphin and Goldman (2013); Bukov and Polkovnikov (2014); Price et al. (2016); Mugel et al. (2017). Following up on these studies, the deflection of an atomic cloud as a response to external force was used to experimentally measure the Chern number in a non-electronic system for the first time Aidelsburger et al. (2015).

While Floquet engineering is a highly flexible and powerful technique, it poses several concerns. One of the main open questions is related to the interplay of driving and interactions which can heat up the system to a featureless, infinite-temperature regime according to general considerations D’Alessio and Rigol (2014); Bukov et al. (2015). In particular, it is shown that an initial Bose-Einstein condensate in a periodically driven optical lattice may become unstable due to two-body collisions Choudhury and Mueller (2015) or through the mechanism of parametric resonance Kennedy et al. (2015); Bukov et al. (2015); Lellouch et al. (2017); Plekhanov et al. (2017); Lellouch and Goldman (2018); Michon et al. (2018); Näger et al. (2018); Boulier et al. (2018). The preparation protocol, stability and a lifetime of strongly correlated phases, expected in the regime of strong interactions under driving is a highly debated topic at the moment Bukov et al. (2015); Lelas et al. (2016); Motruk and Pollmann (2017); Tai et al. (2017).

In order to further explore the role of weak atomic interactions in probing topological features, here we consider the dynamics of weakly interacting incoherent bosons in a driven optical lattice exposed to an external force. The setup that we consider includes all basic ingredients for the Chern-number measurement Dauphin and Goldman (2013); Aidelsburger et al. (2015) - the Chern number of the topological band can be extracted from the center-of-mass motion of atomic cloud in the direction transverse to the applied force. We assume an ideal initial state where the lowest topological band of the effective model is almost uniformly populated. The optimal loading sequence necessary to reach this state is considered in Refs. Ho and Abanin (2016); Dauphin et al. (2017). Following the recent experimental study Aidelsburger et al. (2015), we assume that atoms are suddenly released from the trap and exposed to a uniform force. We perform numerical simulations for the full time-dependent Hamiltonian and take into account the effects of weak repulsive interactions between atoms within the mean-field approximation. We make a comparison between the dynamics governed by the effective and time-dependent Hamiltonian and delineate the contribution of interactions to the center-of-mass response and to the overall cloud expansion dynamics. Our results show that interactions lead to the undesirable atomic transitions between topologcal bands Bilitewski and Cooper (2015), but we also find that a weak atomic repulsion can facilitate the Chern-number measurements in several ways.

The paper is organized as follows: In Sec. II we describe the model and introduce a method that we apply for the description of incoherent bosons. In Sec. III we address the dynamics of noninteracting incoherent bosons, and then in Sec. IV we address the regime of weak repulsive interactions. Finally, we summarize our results in Sec. V. Appendices A to F provide further details.

Ii Model and method

In this Section, we first present the driven model introduced in Ref. Aidelsburger et al. (2015), and then derive the corresponding effective model and discuss its basic characteristics. At the end, we explain our choice of the initial state and outline the method that we use to treat the dynamics of weakly interacting incoherent bosons.

ii.1 Effective Floquet Hamiltonian

Interacting bosons in a two-dimensional optical lattice can be described by the Bose-Hubbard Hamiltonian

(1)

where and are creation and annihilation operators that create and annihilate a particle at the lattice site ( is the lattice constant), is the number operator, and are the hopping amplitudes along and , and is the on-site interaction. In the derivation of the model (1) we use the single-band tight-binding approximation Bloch et al. (2008). Although the experimental setup Aidelsburger et al. (2015) is actually three-dimensional, with an additional confinement in the third direction, our study is simplified to a two-dimensional lattice.

In order to engineer artificial gauge field in the experiment Aidelsburger et al. (2015), hopping along was at first inhibited by an additional staggered potential

(2)

and then restored using resonant laser light. The experimental setup can be described by a time-dependent Hamiltonian

(3)

where is a time-dependent modulation

(4)
(5)

is the driving amplitude and is the resonant driving frequency. We set the relative phase between the optical-lattice potential and the running-waves used for laser-assisted tunneling to .

Figure 1: Schematic representation of the model. The unit cells are shaded. (a) Effective Hamiltonian without correction, (II.1). Vertical links correspond to real hopping amplitudes (along direction), while the horizontal links to the right of lattice sites labeled , , and correspond to complex hopping amplitudes with phases , , and , respectively (when hopping from left to right). (b) Effective Hamiltonian with correction, (II.1). Red lines represent positive next-nearest-neighbor hopping amplitudes (connecting uppercase letters), while the blue lines represent negative next-nearest-neighbor hopping amplitudes (connecting lowercase letters). Nearest-neighbor hopping amplitudes are the same as in (a).

Using Floquet theory, the time-evolution operator corresponding to the Hamiltonian (3) can be represented as

(6)

where is the full time-independent effective Hamiltonian that describes slow motion and is the time-periodic kick-operator that describes micromotion Goldman and Dalibard (2014); Goldman et al. (2015).

For the moment, in this subsection we first consider the non-interacting model . We also assume that the driving frequency is the highest energy scale, but that it is still low enough that the lowest-band approximation used in deriving Eq. (1) is still valid. In the leading order of the high-frequency expansion, the effective Hamiltonian is given by

(7)

where the renormalized hopping amplitudes are and . A schematic representation of this model is presented in Fig. 1(a). The unit cell is shaded and the full lattice is spanned by the vectors and . Particle hopping around a plaquette in the counterclockwise direction acquires a complex phase and the model is equivalent to the Harper-Hofstadter Hamiltonian Hofstadter (1976) for the case Hofstadter (1976). The explicit form of the kick operator from Eq. (3) is given in Appendix A.

Following Refs. Goldman and Dalibard (2014); Goldman et al. (2015), we find that additional corrections of the order contribute to the system’s dynamics and we introduce another approximation for the effective Hamiltonian

(8)

The derivation of Hamiltonian (II.1) is given in Appendix A and its schematic representation is given in Fig. 1(b). The correction introduces next-nearest-neighbor hopping along direction with opposite signs for lattice sites with either even or odd -coordinate . This term does not change the total complex phase per plaquette, but the unit cell is now doubled and thus the first Brillouin zone is halved. A similar term was engineered on purpose in order to implement the Haldane model Jotzu et al. (2014).

In the next subsection we investigate properties of energy bands of both effective Hamiltonians, and . We use the units where and . Unless otherwise stated, we set the parameters to following values: lattice size sites, hopping amplitudes and the driving amplitude . This value of the driving amplitude was chosen to be the same as in the experiment Aidelsburger et al. (2015). In order to set the renormalized hopping amplitude along to , the initial hopping amplitude has to be , and the correction term is therefore proportional to , so it cannot be safely neglected unless the driving frequency is very high.

ii.2 Band structure

Momentum-space representations of the effective Hamiltonians and , denoted by and respectively, are derived in Appendix B. Band structures for the effective Hamiltonian without the correction, Eq. (34), as well as for the effective Hamiltonian including the correction term, Eq. (35), are shown in Fig. 2 for the two values of driving frequencies and .

Figure 2: Energy bands of the effective Hamiltonians. (a) Eq. (34), which is without the correction term. (b) Eq. (35), which includes the correction term. Driving frequency , gaps are open. (c) Same as (b), but with . Gaps are closed.

The Hamiltonian is the Harper-Hofstadter Hamiltonian for the flux . It has four energy bands, where the middle two bands touch at and can therefore be regarded as a single band, see Fig. 2(a). The topological content of these bands is characterized by the topological index called the Chern number. The Chern number is the integral of the Berry curvature Berry (1984) over the first Brillouin zone divided by ,

(9)

where denotes the band number and the Berry curvature is , expressed in terms of eigenstates of the effective Hamiltonian . The Chern numbers of the three well-separated bands are , and .

Because the correction from Eq. (II.1) includes next-nearest-neighbor hopping terms, the elementary cell in real space is doubled (see Fig. 1(b)), and as a consequence, the first Brillouin zone for the Hamiltonian is reduced by a factor of compared to . There are now eight lattice sites in the unit cell and eight energy bands, but the number of gaps depends on the driving frequency. The new bands touch in pairs, in such a way that there are always maximally three well-separated bands. When the driving frequency is high enough, the correction is small and the gaps between the three bands remain open, see Fig. 2(b). The original band structure of is recovered in the limit . The Berry curvature and the Chern number can be calculated using the efficient method presented in Ref. Fukui et al. (2005). Our calculations confirm that the Chern numbers of are equal to those of (, and ), as long as the gaps between the energy bands are open. The gaps close when the driving frequency is too low, see Fig. 2(c), and the Chern numbers of the subbands can no longer be properly defined.

ii.3 Dynamics of incoherent bosons

We need to take into account a contribution of weak, repulsive interactions. Full numerical simulations of an interacting many-body problem are computationally demanding, so we need a reasonable, numerically tractable approximation. To this end we will use the classical field method Kagan and Svistunov (1997) which belongs to a broader class of truncated Wigner approaches Polkovnikov (2010). This method is similar to the approach used to treat incoherent light in instantaneous media Buljan et al. (2004); Cohen et al. (2006), known in optics as the modal theory.

The underlying idea of the method is to represent the initial state as an incoherent mixture of coherent states , , Kagan and Svistunov (1997). This is explained in more details in Appendix C. In our study, we sample initial configurations of these coherent states with

(10)

where are random phases and the states correspond closely to the lowest-band eigenstates of . Each of initial states is time-evolved and physical variables can be extracted by averaging over an ensemble of different initial conditions.

The time evolution of each of these coherent states is governed by

(11)

where are matrix elements of from Eq. (3), is the external force and interactions contribute with the last, nonlinear term. Formally, Eq. (11) takes the form of the Gross-Pitaevskii equation Dalfovo et al. (1999); Pitaevskii and Stringari (2003); Pethick and Smith (2008). The performances and limitations of the method are discussed and reviewed in Ref. S. Gardiner and Szymanska (2013).

For comparison, we also consider the related time evolution governed by the effective Hamiltonian

(12)

where , with being either from Eq. (II.1), or from Eq. (II.1). The equation (12) should be considered only as a tentative description of the system: the mapping between to is strictly valid only in the noninteracting regime and the interaction term may introduce complex, non-local, higher-order corrections D’Alessio and Rigol (2014). However, we expect their contribution to be small in the limit , and for time scales which are not too long Mori et al. (2016); Abanin et al. (2017).

In the following we use modes and accommodate particles per mode, so in total in the simulations we have bosons. Typical densities in real-space are up to particles per site and we choose the values of in the range . Other parameters: , , , . The correction terms are non-negligible in this frequency range. In practice, we first numerically diagonalize the Hamiltonian (37) from Appendix C and set our parameters in such a way that the lowest modes have high overlap with the lowest-band of the effective model. In the next step, we sample initial configurations (10). For each of sets of initial conditions we then time-evolve Eq. (11) and extract quantities of interest by averaging over resulting trajectories. This value of is chosen to be high enough, so that the fluctuations are weak. We present and discuss results of our numerical simulations in the following sections.

Iii Noninteracting case

We start by addressing the dynamics of noninteracting bosons. In this case we set in Eq. (11) and numerically solve the single-particle Schrödinger equation without further approximations. Our aim is to numerically validate and compare the two approximations, Eqs. (II.1) and (II.1) for the effective Hamiltonian. To this purpose, we juxtapose results of the two approximative schemes with the numerically exact results obtained by considering the full time evolution governed by . For clarity, the four different time evolutions that we consider in this section are summarized in Table 1. We calculate the center-of-mass position and plot the results in Fig. 3. In this way we also find the regime of microscopic parameters where the Chern-number measurement can be optimally performed.

First, we consider the basic Harper-Hofstadter Hamiltonian (II.1) and select the occupied modes of the initial state (36) as eigenstates of the model from Eq. (10) for . As explained in the previous section, at the initial moment , the confinement is turned off and the force is turned on. As a consequence of the applied external force and the non-zero Chern number of the lowest band of the model (II.1), the particles exhibit an anomalous velocity in the direction perpendicular to the force Xiao et al. (2010). In the ideal case, when the lowest band is fully populated, the theoretical prediction for the center-of-mass position in the direction is Aidelsburger et al. (2015)

(13)

where is the Chern number (9) of the lowest band. However, even in the ideal case, due to the sudden quench of the linear potential, a fraction of particles is transferred to the higher bands. To take this effect into account, the authors of Ref. Aidelsburger et al. (2015) introduced a filling factor

(14)

where are populations of different bands of Hamiltionian (II.1) from Eq. (39) in Appendix C and the plus and minus signs in Eq. (14) are defined according to the Chern numbers , and . The final theoretical prediction is then Aidelsburger et al. (2015)

(15)
case initial state band populations evolution
1
2
3
4
Table 1: Four different cases: the same effective Hamiltonian is always used for the initial state and band definitions/calculation of band populations, either with or without the correction. The evolution is governed either by the time-dependent Hamiltonian or by the same effective Hamiltonian as the one that was used for the initial state and band populations.
Figure 3: Anomalous drift . Dashed purple lines: numerical simulations using the time-dependent Hamiltonian (cases 2 and 4 from Table 1). Solid green lines: effective Hamiltonians (c) and (d), and (a) and (b) (cases 1 and 3). Dotted black lines: theoretical prediction (15) from or . (a) Initial states and band populations obtained using the effective Hamiltonian without the correction (cases 3 and 4). Driving frequency . (b) . (c) Hamiltonian with the correction (cases 1 and 2). Driving frequency . (d) .

In Fig. 3(a) we consider the anomalous drift for a high value of the driving frequency , where we expect the expansion in to be reliable. We find an excellent agreement between the prediction (15) (dotted black line) and numerical calculation based on (solid green line). However, some deviations between the full numerical results (dashed purple line) and the results of the approximation scheme (solid green line) are clearly visible. These deviations are even more pronounced for , Fig. 3(b).

Now we turn to the effective model (II.1). In this case we select the modes of the initial state as eigenstates of Eq. (10) for . Moreover, we also consider band populations (39) of the same model. In the case when , Fig. 3(c), the anomalous drift obtained using the effective Hamiltonian (II.1) (solid green line) closely follows the theoretical prediction (15). Moreover, from the same figure we can see that the effective Hamiltonian reproduces the behavior of the time-dependent Hamiltonian very well. All three curves almost overlap for intermediate times (), see Fig. 3(c). We attribute the long-time () deviations to the finite-size effects introduced by the next-nearest-neighbor hopping terms, which cause the atomic cloud to reach the edge of the lattice faster. This effect is explained in more detail in section IV.2.

For a lower driving frequency , the effective and the time-dependent Hamiltonians do not agree so well anymore, see Fig. 3(d). The finite-size effects can be observed even earlier in this case (around ). This happens because the next-nearest-hopping terms are inversely proportional to the driving frequency. It is interesting to note that the prediction (15) is close to numerical data for short times even in this case when the gaps of the effective model are closed, see Fig. 2(c), and the Chern number of the lowest band is not well defined. In fact, it is surprising that the anomalous drift even exists in this case, as all subbands are now merged into a single band.

Figure 4: Time-evolution of the filling factor for driving frequency . Solid purple lines: evolution governed by the time-dependent Hamiltonian (cases 2 and 4 from Table 1). Dashed green lines: evolution governed by the effective Hamiltonian or (cases 1 and 3). Dotted black lines: green lines shifted in order to compare them with purple lines. Shift is chosen so that the two lines approximately overlap. (a) Initial states and band populations obtained using the effective Hamiltonian , which is without the correction term (cases 3 and 4). (b) Hamiltonian which is with the correction term (cases 1 and 2).

Time-evolution of the filling factor is plotted in Fig. 4 for four different cases from Table 1 – evolution using the effective Hamiltonian without correction (, case 3, dashed green line in Fig. 4(a)) the effective Hamiltonian with correction (, case 1, dashed green line in Fig. 4(b)), or the time-dependent Hamiltonian (, cases 2 and 4, solid purple lines). At the initial moment , because the initial state was multiplied by the operator . This introduces a shift between and . Apart from the shift, these two curves behave similarly, unlike the curve that exhibits completely different behavior. Because of this, we use only to estimate the value of the prediction (15).

We find that the values of for are high , Fig. 4. For this reason, up to the center-of-mass position exhibits roughly linear behavior with some additional oscillations. Interestingly, the anomalous drift exhibits quadratic behavior on short time scales in all cases from Fig. 3. In Appendix D, we explain this feature using the time-dependent perturbation theory and the Fermi’s golden rule.

Iv Interacting case

We now investigate the effects of weak repulsive interactions. We work in the high frequency regime and set . As shown in subsection II.2, for the effective Hamiltonian with correction, , is in this case equivalent to the Harper-Hofstadter Hamiltonian with flux . Moreover, the same approximative form of the full effective model accurately reproduces the behavior of the time-dependent Hamiltonian up to and thus provides a good starting point for the study of weakly interacting particles. We first consider the anomalous drift of the center-of-mass of the atomic cloud and then we inspect the expansion dynamics more closely in terms of atomic density distributions in real and momentum space.

iv.1 Anomalous drift and dynamics of band populations

To simulate the dynamics of many incoherent bosons, we use the classical field method presented in subsection II.3 and propagate Eq. (11) in time. We assume that at atoms are uniformly distributed over the lowest band of . For this reason, the initial state is the same as the one that we use in the noninteracting regime. In this way, the dynamics is initiated by an effective triple quench: at the confining potential is turned off, atoms are exposed to the force and also the interactions between particles are introduced. The total number of particles is set to , which amounts to approximately particles per lattice site in the central region of the atomic cloud. We consider only weak repulsion .

The anomalous drift obtained using the full time-dependent Hamiltonian is shown in Fig. 5(a) for several different values of the interaction strength . In comparison to the noninteracting regime, we find that the weak repulsive interactions inhibit the response of the center-of-mass to the external force. In particular, at the drift is reduced by about for and it is further lowered by an increase in . Finally, at , the anomalous drift is barely discernible. Interestingly, for weak we find that the drift in the range of looks “more linear” as a function of time in comparison to the noninteracting result.

Figure 5: The effects of interactions. (a) Anomalous drift for several different values of the interaction coefficient . is given in units where . Full lines – numerical simulations using the time-dependent Hamiltonian . Dashed lines – theoretical prediction (15) from . (b) Corresponding , obtained from simulations using the effective Hamiltonian .

We now analyze the anomalous drift in terms of the filling factor and compare the results of Eq. (11) with the description based on Eq. (12). By solving Eq. (12) we obtain the filling factor following Eq. (39) and present our results in Fig. 5(b). Whenever the results of Eq. (11) reasonably agree with the results obtained from Eq. (12), we are close to a steady-state regime with only small fluctuations in the total energy, as Eq. (12) preserves the total energy of the system. In this regime, during the expansion dynamics the interaction energy is converted into the kinetic energy and atoms are transferred to higher bands of the effective model. Consequently, the filling factor is reduced. Typically, we find three different stages in the decrease of .

In an early stage, , a fast redistribution of particles over the bands of the effective model sets in due to the sudden quench of . The factor decays quadratically as a function of time down to for , and for . In this process the interaction energy of the system is quickly lowered as described in Appendix E. At later times , we observe a linear decay of the filling factor as a function of time, that finally turns into an exponential decay at even later times (). Similar regimes are observed in other dynamical systems. For example, a decay rate of an initial state suddenly coupled to a bath of additional degrees of freedom exhibits these three stages Debierre et al. (2015). The initial quadratic decay is often denoted as “the Zeno regime”. For longer propagation times, the Fermi’s golden rule predicts the linear decay. At even longer time scales, when the repopulation of the initial state is taken into account, the time-dependent perturbation theory yields the exponential regime, known under the name of the Wigner-Weisskopf theory Debierre et al. (2015).

We now investigate this last regime in more detail. For the population of the lowest band , an exponential decay function provides high quality fits for , see Fig. 6(a) for an example. Similarly, the populations of two higher bands can also be fitted to exponential functions. The obtained exponential decay coefficients for the lowest band population are plotted as a function of the interaction strength in Fig. 6(b). The resulting dependence is approximately quadratic: . For small values of , the exponents obtained for the dynamics governed by and agree very well and exhibit linear behavior. At stronger interaction strengths , the approximation of Eq. (12) becomes less accurate as it omits the quadratic contribution in found in the full time evolution. In addition, the values of the exponents are affected by the force strength and driving frequency .

Figure 6: (a) Effective Hamiltonian evolution of the band populations . Full lines – numerical results. Dashed lines – exponential fit using . The coefficient was fixed to , and for the first, second and third band respectively. (b) Dependence of the exponential decay coefficients for the lowest band population on the interaction strength. is given in units where .

As we now understand some basic features of , we make an explicit comparison between the numerical results for the anomalous drift and the expectation (15). The dashed lines in Fig. 5(a) correspond to the theoretical prediction (15) calculated from . For the intermediate interaction strengths , we find a very good agreement between the two. From this we conclude that the interaction-induced transitions of atoms to higher bands are the main cause of the reduced anomalous drift as a function of . When the interactions become strong enough (), the numerical results start to deviate from the theoretical prediction (15) with . In this regime, Eq. (12) does not provide a reliable description of the dynamics, as higher order corrections need to be taken into account.

iv.2 Real and momentum-space dynamics

Figure 7: Real-space density distribution, noninteracting case . (a) Initial state. (b) After ( driving periods), evolution using the time-dependent Hamiltonian . (c) Evolution using effective Hamiltonian without correction . (d) Evolution using effective Hamiltonian with correction .

At the initial moment, the atomic cloud is localized in the center of the lattice. By setting in the confining potential of Eq. (37) and populating the lowest-lying states, we fix the cloud radius to , Fig. 7(a). The cloud density is of the order of atoms per lattice site and a weak density modulation is visible along direction. After the confining potential is turned off, and the external force in the direction is turned on, the cloud starts to expand and move in the direction. As shown in the previous subsection, the band populations and therefore the anomalous drift are significantly altered by the interaction strength, and this is also the case with the expansion dynamics, see Figs. 7 and 8.

In the noninteracting case, Fig. 7(b), the atomic cloud nearly separates into two parts moving in opposite directions along axes (while the center of mass still moves in the direction). By comparing Fig. 7(c) and Fig. 7(d), we conclude that this effect stems from the next-nearest-neighbor hopping along present in the effective Hamiltonian (II.1), as it does not happen in the effective model without the correction term (II.1). This type of separation was already observed in Ref. Dauphin and Goldman (2013), where the next-nearest-neighbor hopping terms were also present.

Figure 8: Real-space density distribution after ( driving periods), interacting case. is given in units where . (a) Evolution using the time-dependent Hamiltonian , . (b) Same with . (c) Evolution using the effective Hamiltonian , . (d) Same with .

So far we have considered the averaged response of the whole atomic cloud. We now inspect the expansion dynamics in a spatially resolved manner. The real-space probability densities at the initial moment and after ( driving periods) are shown in Figs. 7 and 8, and the corresponding momentum-space probability densities in Appendix F.

Figure 9: Atomic cloud width for different interaction strengths, evolution using the time-dependent Hamiltonian . is given in units where . (a) . (b) .
Figure 10: (a) Comparison of anomalous drifts obtained from evolution using the time-dependent Hamiltonian (solid purple line), effective Hamiltonian without correction (dashed green line) and effective Hamiltonian with correction (dotted black line). Intermediate interaction strength . is given in units where . (b) Time-evolution of the inverse participation ratio in momentum space for several different values of . Evolution is performed using the time-dependent Hamiltonian . When the interactions are strong enough, approaches the maximal possible value ( in this case), which is equal to the total number of states and corresponds to the completely delocalized state. is given in units where . (c) Chern number of the lowest band obtained for different interaction strengths as the ratio of the theoretical prediction for the anomalous drift and numerical results: .

When the interactions between particles are included, this separation is not so prominent (Fig. 8(a), ), and it almost completely disappears when the interactions are strong enough (Fig. 8(b), ). This is also the case when the evolution is governed by the effective Hamiltonian , see Figs. 8(c) and 8(d). Atomic cloud widths during the expansion are plotted in Fig. 9. We observe a slow expansion of the cloud in direction, Fig. 9(b), and much faster expansion along direction, Fig. 9(a), which comes about as a consequence of the cloud separation. On top of this, we observe that the interactions enhance expansion along . Surprisingly, the opposite is true for the dynamics along . This counterintuitive effect is often labeled as self-trapping and its basic realization is known for the double-well potential Milburn et al. (1997); Raghavan et al. (1999). In brief, strong repulsive interactions can preserve the density imbalance between the two wells, as the system can not release an excess of the interaction energy. In our case, the situation is slightly more involved as the cloud splitting is inherent (induced by the corrections of the ideal effective Hamiltonian). Apart from this, due to the driving the total energy is not conserved. However, our numerical results indicate that the interaction energy is slowly released in the second expansion stage, Fig. 13. Effectively, in this way the interactions cancel out the contribution of the next-nearest neighbor hopping and favor the measurement of the properties of the model (II.1). In Fig. 10(a) we show that deviations between different approximations based on , and in the anomalous drift nearly vanish at .

Another desirable effect might be that the interactions make the momentum-space probability density more homogeneous, see Appendix F, so that the real-space probability density becomes more localized. We can quantify momentum-space homogeneity using the inverse participation ratio , where is the probability that the state is occupied at time . Minimal value of is and it corresponds to a completely localized state, while the maximal value is equal to the total number of states (in our case ) and corresponds to the completely delocalized state, where the particles have the same probability of being at any quasi-momentum . As stated before, the first Brillouin zone of the lowest band has to be as homogeneously populated as possible in order to properly measure the lowest band Chern number. From Fig. 10(b), we see that increases in time when the interaction coefficient is large enough, so we can conclude that the interactions are actually beneficial for measuring the Chern number, as they can “smooth-out” the momentum-space probability density. In Fig. 10(c) we give estimates for the Chern number that can be extracted from our numerical data for different values of . We find the best estimate for the intermediate interaction strength .

V Conclusions

Motivated by the recent experimental results reporting the Chern numbers of topological bands in cold-atom setups, we studied numerically bosonic transport in a driven optical lattice. The considered driving scheme and the range of microscopic parameters were chosen to be close to those in a recent experimental study Aidelsburger et al. (2015). The driving frequency was set to be high enough in order to avoid strong energy absorption for the relevant time scales. We investigated bosonic dynamics for the full time-dependent Hamiltonian, the effective Floquet Hamiltonian, and included the effects of weak repulsive interactions between atoms using the mean-field approximation. In the noninteracting case, we found that the effective Hamiltonian and its band structure depend on the frequency of the drive through an additional correction term.

The main focus of this work is on the effects of weak interactions. For a weak atomic repulsion, atomic transitions to higher effective bands obtained in our simulations mainly occur due to a release of the initial interaction energy during the atomic-cloud expansion. Although the effect is undesirable, it can be properly taken into account in the extraction of the Chern number. At larger interaction strengths, the transitions are more pronounced as the system absorbs energy from the drive. In this regime the good agreement between the full and effective description is lost and the measurement should become more complicated. In addition to causing redistribution of atoms over bands, our results show that weak interactions can also be beneficial in measuring the Chern number. Their desirable effect comes about due to smoothening the atomic distribution over the topological band and due to canceling out the contribution of some less relevant terms to the bosonic dynamics.

Acknowledgements.
This work was supported by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under Projects ON171017, BKMH and TOP-FOP, the Croatian Science Foundation under Grant No. IP-2016-06-5885 SynthMagIA and the TOP-FOP project, DAAD (German Academic and Exchange Service) under the BKMH project and by the Deutsche Forschungsgemeinschaft via DFG FOR 2414 and the high-performance computing center LOEWE-CSC. Numerical simulations were performed on the PARADOX supercomputing facility at the Scientific Computing Laboratory of the Institute of Physics Belgrade.

Appendix A The effective model

After a unitary transformation into the rotating frame , where and are the old and the new wave-function, and is the staggered potential, the new time-dependent Hamiltonian that describes the experimental setup is given by Aidelsburger et al. (2015)

(16)

where

(17)
(18)
(19)

The kick-operator is given by

(20)

and the effective Hamiltonian by

(21)

If we assume that the driving frequency is high and interactions are weak, the interaction term and almost all terms can be neglected. After substituting Eqs. (16), (17) and (18) into Eq. (A) we obtain:

(22)
(23)

We will now separately calculate each term:

(24)
(25)
(26)
(27)

Using trigonometric identities and

(28)
(29)

we can rewrite the sum of terms (A) and (A) in a more convenient form

(30)

The only () term that cannot be neglected in the parameter range that we use is Aidelsburger et al. (2015)

(31)

Finally, the effective Hamiltonian becomes

(32)