Topological Phase Diagram of Optimally Shaken Honeycomb Lattices: A Dual Perspective from Stroboscopic and Non-Stroboscopic Floquet Hamiltonians

We present a direct comparison between the stroboscopic and non-stroboscopic effective approaches for ultracold atoms in shaken honeycomb lattices, focusing specifically on the optimal driving introduced by A. Verdeny and F. Mintert [Phys. Rev. A 92, 063615 (2015)]. In the fast-driving regime, we compare the effective non-stroboscopic Hamiltonian derived through a perturbative expansion with a non-perturbative calculation of the stroboscopic Floquet Hamiltonian, obtained through a simple non-perturbative numerical approach. We show that while some of the tunneling parameters are inherently model-dependent, the topological properties of the system remains robust, as expected. Using the same numerical approach we compute the topological phase diagram, arguing that it is most effectively represented in terms of the physical parameters characterizing the driving and the bare Hamiltonian -- parameters directly accessible in experiments -- rather than the emergent tunneling parameters, that depend on the model representation.


I. INTRODUCTION
Recent years have witnessed the flourishing of the socalled Floquet engineering in the context of of ultracold quantum gases, namely the ability to experimentally manipulate and control the behavior of these systems by applying time-periodic driving forces .This approach takes advantage of the Floquet theory, which describes the behavior of quantum systems under the influence of time-periodic perturbations (for a review see, e.g., Refs.[15][16][17]).The key idea behind all this is to engineer the periodic driving force in such a way that it effectively modifies the behavior of the system on timescales much larger than the modulation period.By controlling the parameters of the driving force, one can tailor the energy spectrum of the system, induce artificial gauge fields, create topological states of matter, and realize novel quantum phases that are not possible in static systems.Nowadays, many research laboratories worldwide have the capability to experimentally implement such engineering techniques [2,9,12,25,29,31,38].
From a theoretical perspective, different approaches are available to describe such periodically driven systems.A natural choice is provided by the original Floquet theory, in which the evolution operator is factorized in the product of a time-periodic evolution operator and of a one-cycle evolution operator, see, e.g., Ref. [15].Through the latter one can introduce an effective Floquet Hamiltonian ĤF [t 0 ] [12,24], providing an effective description of the longtime dynamics of the system.It corresponds to the full evolution of the system on a whole period T , from an initial time t 0 to t 0 + T .For this reason, it is also called stroboscopic Floquet Hamiltonian [16].Though ĤF [t 0 ] depends on the choice of the time t 0 which defines the beginning of the stroboscopic driving period, it is otherwise time-independent, in the sense that it does not describe the micromotion within each period T .It is noteworthy that the energy spectrum and topological properties of the effective Hamiltonian do not depend on the choice of the initial time t 0 , since two effective Hamiltonians for different initial times are related through a unitary gauge transformation, ĤF [t ′ 0 ] = Û (t ′ 0 , t 0 ) ĤF [t 0 ] Û −1 (t ′ 0 , t 0 ).The choice of the initial time t 0 only affects the local properties of the Hamiltonian, i.e., the actual values of the effective tunneling coefficients, but not the general structure of the Hamiltonian.In principle, since the choice of t 0 is a gauge choice, all the Hamiltonians ĤF [t 0 ] must be gauge equivalent to some fixed Floquet Hamiltonian ĤF , that is independent of t 0 [16].
In general, an exact expression for the effective Hamiltonian ĤF [t 0 ] cannot be found.However, from its formal definition, ĤF [t 0 ] ≡ (i/T ) ln[ Û (t 0 + T, t 0 )], it can be calculated either perturbatively or numerically.The perturbative approach amounts to the so-called Magnus expansion (see, e.g., Refs.[16,39]), that is an analytic Taylor expansion of ĤF [t 0 ] in powers of ω −1 (with ω ≡ 2π/T ).On the numerical side, a convenient approach -that will be employed in this paper -consists in factorizing the time-ordered product over infinitesimal time intervals and handling the resulting expression in the basis of the Pauli matrices [12,24].
A different but related strategy is the one introduced by Goldman and Dalibard in Ref. [13] (see also [16] and references therein), in which the evolution operator is factorized in three terms, one corresponding to the effects associated with the initial phase of the modulation, another one to the longtime dynamics of the systemto which one can associate an effective Hamiltonian Ĥeff -, and finally one accounting for the micromotion.No-tably, Ĥeff is independent of t 0 by definition.Therefore, with respect to the standard Floquet approach described above, it offers the advantage that all quantities are manifestly gauge invariant, by definition.However, in general it cannot be defined by an explicit analytical expression but only in terms of a perturbative expansion in powers of the inverse driving frequency ω −1 [13].The latter is the analogue of the Magnus expansion for ĤF [t 0 ].Indeed, one can establish a formal correspondence between ĤF [t 0 ] and Ĥeff , at each order of the perturbative expansion in ω −1 , see, e.g., Ref. [16].The other side of the coin is that without an explicit analytical expression defining the Hamiltonian Ĥeff , it is generally not possible to compute it numerically in a non-perturbative manner.
To gain a deeper understanding of these concepts, in this paper we provide a direct comparison between the two approaches, by considering a specific illustrative example.To this end, we chose the proposal by Verdeny and Mintert for realizing a tunable Chern insulator [20], consisting in a system of noninteracting ultracold atoms loaded in a honeycomb lattice and subjected to a bichromatic optimal driving [40].As in the original proposal, the model is considered here within the tight-binding approximation.This system presents interesting similarities with the Haldane model [41], a paradigmatic lattice model of a topological insulator, whose phase diagram has been recently explored in experiments with ultracold fermionic atoms in periodically driven optical lattices [12] (see also Ref. [24]).
Specifically, here we focus the analysis on the fastdriving regime, in which the first-order expression of Ĥeff derived in Ref. [20] is expected to be fully justified.We compute the effective (stroboscopic) Floquet Hamiltonian ĤF [t 0 ] by means of the non-perturbative numerical approach of Ref. [24], that allows for a straightforward calculation of the tunneling coefficients and of the topological phase diagram.These results are then compared to the first-order prediction for Ĥeff [20].In particular, we show that within the parameter regime considered in this paper, the perturbative approach well reproduces the topological properties of the model.We also discuss how to conveniently draw the phase diagram in terms of the model parameters associated with the driving and the explicit breaking of the inversion symmetry (IS), without relying on the specific form of the first-order effective Hamiltonian.This is done while considering the fact that IS is dynamically broken by this class of optimal driving.Our discussion also encompasses other interesting aspects, such as the gauge and model dependence of the tunneling coefficients, which, however, is not relevant to the topological properties.
The paper is organized as follows.In Sec.II, we outline the theoretical framework for driven lattice models, specifically focusing on the honeycomb lattice with nearest-neighbor hoppings under the influence of bichromatic sinusoidal driving.In Sec.III, we discuss the stroboscopic and non-stroboscopic approaches for computing the effective Floquet Hamiltonian.First, in Sec.III A, we revisit the perturbative formulation for constructing the 'optimal' effective Hamiltonian as presented in Ref. [20].Then, in Sec.III B, we present the methodology for constructing the stroboscopic Floquet Hamiltonian using a non-perturbative numerical approach.We also provide results for the onsite energies and tunneling coefficients of the model, along with a discussion of analogies and differences between the two approaches.Section IV is dedicated to the discussion of the topological phase diagram.Finally, we summarize our findings and present concluding remarks in Sec.V.

II. DRIVEN LATTICE MODELS
Let us consider a general undriven lattice model described by the bare Hamiltonian where the indices j, j ′ run over the unit cells of the Bravais lattice while the indices ν, ν ′ run over the different sites of a given cell.The coefficients t j−j ′ νν ′ represent the hopping matrix elements in the model, namely the tunneling coefficients.
If the lattice is subjected to a rigid time-periodic driving r lat (t), a particle of mass m loaded in the lattice feel an inertial force F (t) = −mr lat ≡ q(t), with F (t + T ) = F (t). Therefore, the total Hamiltonian of the system may be expressed as [12,24] that, in the co-moving frame, can be conveniently written in the compact form as where and under the condition of vanishing net transferred momentum (over a period), with t b being an arbitrary integration bound.In the following we shall work in this representation (namely, the co-moving frame), omitting the index lat for easiness of notations.We can also fix t b = 0, without loss of generality.
When one is interested in the dynamics of the system over a time scale much larger than the driving period T , an effective time-independent Hamiltonian may be introduced in terms of the one-period time-evolution operator, the so called effective Floquet Hamiltonian (see, e.g., Ref. [24] and references therein), where T represents the time-ordering operator.The above expression should be interpreted in the context of its Taylor expansion.Specifically, by introducing the dimensionless time parameter τ ≡ ωt, the evolution operator can be elegantly expressed as follows, representing a series expansion in 1/ω, commonly known in the literature as the Magnus expansion [16,42].Notice that the operator ĤF [t 0 ] is invariant under the transformation Û (t 0 , t 0 +T ) −→ Û (t 0 , t 0 +T )e i2mπ (with m ∈ Z), so that its eigenvalues, the so called quasienergies, are defined modulo an integer multiple of ω [15].We also recall that the energy spectrum and topological properties do not depend on the initial time t 0 , which only affects the local properties of the Hamiltonian.Owing to this, we can set t 0 = 0, without loss of generality.The dependence on t 0 of the tunneling coefficients within the stroboscopic approach will be briefly discussed later on, in Sec.III B.
In the following, we shall focus on the model considered in Ref. [20], that we review below.

A. Optimally shaken honeycomb lattices
Bare Hamiltonian.The undriven system, represented by the Hamiltonian Ĥb [see Eq. ( 1)], consists of a honeycomb lattice with two inequivalent sites per unit cell as shown in Fig. 1.The geometry of the lattice is identified by two sets of vectors, a i and b i (i = 1, 2, 3), that connect nearest-neighbour (NN) sites (of type A and B) and neighbouring unit cells, respectively.The model is characterized by an onsite energy offset 2δ (where the onsite energies are E A,B = ±δ) and an isotropic NN hopping matrix element between sites A and B, represented by the tunneling coefficient j 0 .All higher order tunneling coefficients are assumed to be vanishing (they are supposed to be negligible, from the physical point of view).
In reciprocal space (see Appendix A), the system Hamiltonian is represented by a 2 × 2 matrix h νν ′ (k) that can be conveniently written in a compact form, by using the basis formed by the 2 × 2 identity matrix, I, and of the three Pauli matrices, σ i [20,24].Namely, with h ≡ (h 1 , h 2 , h 3 ).In the present case, owing to the vanishing of the tunneling coefficients beyond NN and to the fact that E A + E B = 0, we have h 0 (k) ≡ 0. Remarkably, we shall see that this property -that constitutes a distinctive difference with respect to the original Haldane model -is preserved in the presence of the driving, so that h 0 (k) can be omitted from the following discussion.Therefore, the spectrum of h(k) can be readily written as In reciprocal space, the bare Hamiltonian has the following structure where the function is fixed by the lattice geometry.Lattice driving.To induce topologically nontrivial phases, a time-reversal symmetry (TRS) breaking driving must be applied.Here, we specifically consider the bichromatic driving proposed in Ref. [20], which is designed to achieve isotropic effective nearest-neighbor Trajectory of the lattice (r lat = −F (t)/m) under the influence of the force (14), for some values of the ratio A2/A1.The path is given over a whole period; the thick dot at (x, y) = (0, 0) indicates the initial position, at t = 0. Lengths (x, y) are given in units of A1/(mω 2 ).The orientation of the axis is the same as in Fig. 1.Notice that when A2 = 0, IS (with respect to the trajectory center) is always explicitly broken.However, the trajectory profiles always exhibit at least C3v symmetry.
(NN) and next-nearest-neighbor (NNN) tunneling rates.Namely, where e x,y are the unit vectors of the x and y axis and ω = 2π/T is the fundamental frequency of the driving.Typical trajectories of the lattice center-of-mass, produced by such a force, are shown in Fig. 2. Note that this driving scheme differs from the monochromatic driving employed by Jotzu et al. [12] in their experimental realization of the Haldane topological phase diagram.In that case, the breaking of TRS was directly controlled by the phase ϕ of the driving (see also Ref. [24]).
As anticipated in the introduction, it is important to remark that this class of driving explicitly breaks IS.This is the reason why the model can exhibit a Haldanelike topological phase diagram, even if the undriven lattice Hamiltonian is inversion-symmetric (see Ref. [20]).Bearing this in mind, in the following we will refer to the inversion-symmetric (δ = 0) and IS-breaking (δ = 0) cases in reference to the properties of the undriven model, specifically whether the bare onsite energies of sites A and B are degenerate or not.
From Eq. ( 4) it becomes apparent that the driving leaves the onsite energies unaffected, t 0 AA (t) = δ = −t 0 BB (t), whereas the NN hopping elements t j−j ′ νν ′ (t) acquire a time dependence.In particular, by defining [20] g ai (t) ≡ j 0 e iχa i (t) , with the three NN tunneling coefficients corresponding to the hopping from the site A in the central unit cell to three neighbouring sites B (see Fig. 1) can be expressed as By combining the above results, the components of the driven Hamiltonian in reciprocal space read with h 0 (k, t) ≡ 0.

III. EFFECTIVE FLOQUET HAMILTONIAN
In general, an exact expression for the effective Hamiltonian, ĤF [t 0 ] in Eq. ( 6) or Ĥeff [13,16], is unavailable.However, it can be calculated either perturbatively or numerically, as previously mentioned.Different analytical approaches are available in the literature, as discussed in Ref. [16].However, perturbative expansions are often limited to the first order and this may not be sufficient for an accurate quantitative description in a generic driving regime.In this respect, a non-perturbative numerical approach can be very useful.In the following, we will directly compare these two viewpoints within the framework of the driven lattice model presented in the previous section.

A. Perturbative approach
We begin by examining the perturbative approach, which amounts to an expansion of Ĥeff in powers of ω −1 .In particular, we adopt the approach of Verdeny and Mintert [20], which serves as a formal definition of the effective Hamiltonian, where the operators Ĥn are the Fourier components of the time-dependent Hamiltonian (3), This expression is obtained within the approach discussed by Goldman and Dalibard [13], differing from the usual Magnus expansion by lacking additional commutator terms (see, e.g., Ref. [16]).We recall that since Ĥ(0) eff = Ĥ0 is defined as the time average of Ĥ(t), it is characterized by the same physical processes of the undriven Hamiltonian, whose coefficients simply get renormalized, as we shall discuss in the following discussion.
Instead, the structure of Ĥ( 1) eff (which involves a commutator) leads -already at first order -to the emergence of new couplings (that is, new tunneling processes).This, in essence, is the reason why "Floquet engineering" is such a powerful tool for realizing non trivial topological models.
Let us then consider the effective Floquet Hamiltonian defined by Eq. ( 6), starting with its first-order truncation in Eq. ( 21).As discussed in Ref. [20], a bichromatic driving designed as in Eq. ( 14) yields isotropic NN and next-nearest-neighbours (NNN) tunnelling rates, within the first order approximation in Eq. ( 21).The leading order term Ĥ(0) eff = Ĥ0 is obtained from Eq. ( 22).Therefore, it is straightforward to see that the onsite energies remain unchanged (they do not depend on time), whereas the NN hopping elements get renormalized by time-averaging the expressions in Eq. ( 17), e.g.(we omit the label 'eff' for easiness of notations), As for the first order term Ĥ(1) eff , it is characterized by both a renormalization of the onsite energies and, more interestingly, by the emergence of NNN hopping terms, which were vanishing in the undriven model.In the above expressions, the indexes i, j, k are cyclic integers, for example t b2 AA = β(a 3 , −a 1 ), and the function β(a i , a j ) is defined as follows, where the functions g n aj are the Fourier components of Eq. ( 15) , i.e., g n aj = T −1 T 0 dt g aj (t)e −inωt [20].Intuitively, one may interpret the function β as the rate for a particle to tunnel to a site by the vector a i and then to tunnel again through the vector a j [19].As anticipated above, it can be demonstrated [20] that the NNN coefficients turn out to be isotropic, namely they are equal along the tree directions.Therefore, the NN and NNN tunneling rates can be rewritten in a compact form, with the superscript in parenthesis denoting the order of the inverse frequency expansion the parameters are related to, as (1) Then, the effective Hamiltonian h eff (k) reads The above results, obtained by Verdeny and Mintert in Ref. [20], will serve as a reference for comparing the results of the numerical approach discussed in the following section.As an example, in Fig. 3 we show the behavior of the three t.b. parameters, ∆, t (0) 0 , and t (modulus and phase [43]), as a function of the driving amplitudes.In particular, we consider the case in which there is no explicit IS breaking, δ = 0, under fast driving conditions, with ω = 100j 0 / .This high-frequency value is expected to be sufficiently large to warrant the first-order truncation of the inverse frequency expansion in Eq. ( 21).

B. Non-perturbative numerical approach
In this section, we will extend the analysis beyond the lowest order in the inverse frequency expansion.To this end, we will employ the numerical approach discussed in Ref. [24].This approach corresponds to summing up the entire inverse frequency expansion, allowing us to compute the effective Floquet Hamiltonian ĤF without additional approximations other than those required by the numerical implementation, such as time discretization.We recall that ĤF is expected to be gauge-equivalent to Ĥeff [16].This implies that the two Hamiltonians should share the same spectrum and topological prop- erties.However, 'local' properties, such as the tunneling coefficients, may not necessarily be equal.
To start with, we consider the effective Floquet Hamiltonian in reciprocal space, h F (k), that we compute by means of a suitable time discretization, as follows with ∆ t ≡ T /N .The details of numerical scheme are summarized in Appendix B. Then, the tunneling coefficients can be easily obtained as see Appendix A. Consistently with the notations in Sec.III A, we indicate the non-perturbative (np) effective onsite energies, as well as the NN and NNN tunnelling rates obtained through this method as ∆ np ν , t np 0i , and t np 1ν,i = |t np 1ν,i |e iφ np νi (ν = A, B), respectively.The index i, introduced to account for any possible non-isotropy, matches that of the vectors a i in case of NN tunneling rates (i = 1, 2, 3), and of b i for NNN hoppings (with i = 4, 5, 6 for −b i ), see Fig. 1.
The behavior of the various tight-binding parameters of the model as a function of the amplitudes A 1 and A 2 of the driving is discussed below.We focus on the same scenario as considered in Section III A, specifically addressing the case where there is no explicit IS breaking, δ = 0, under fast driving conditions, with ω = 100j 0 / .
We begin by illustrating the behavior of the onsite energies and of the NN tunneling coefficients.Regarding the former, we observe that the two lattice sites exhibit opposite onsite energies, as obtained in the model of Ref. [20], see Eq. (24).Therefore, moving forward, we can simplify our notation by setting ∆ np ≡ ∆ np A = −∆ np B .This quantity is shown in Fig. 4a, which nicely matches with Fig. 3a.These figures clearly reveal the dynamical breaking of the degeneracy between the two lattice sites caused the driving.
Similarly, we find that the NN tunnelling rate t np 0i is isotropic, as in the non-stroboscopic model, so that the index i can be safely omitted.Its behavior, shown in Fig. 4b, perfectly reproduces that displayed in Fig. 3b, proving that also the emergent NN tunneling is modelindependent.This was expected because in the fastdriving regime the value of the NN tunneling coefficient is determined by the zeroth order effective Hamiltonian, where both the stroboscopic and non-stroboscopic have the same expression as in Eq. ( 22) (with n = 0), and therefore they coincide.This provides a robust characterization of the system's topological properties, as from the isotropy of the NN hoppings it follows that the position of the Dirac points is fixed by the space group of the model, see Appendix C.This has to be a gauge invariant property, so that it cannot depend on the model, as confirmed by the numerical results.Then, in Fig. 5 we show the modulus of the emergent NNN tunnelling rates, |t np 1ν,i |, and their complex phases φ np νi , in the left and right columns respectively.Consistently with the property in Eq. ( 25) found from the perturbative approach [44], the numerical results obey the relation t np 1B,1 = −t np 1A,1 .Therefore, in the figure and in the following discussion the ν index is omitted, for the sake of notation simplicity.
On the other hand, we find that the NNN tunneling rates are model-dependent, both in modulus and phase, see Figs. 3c,d for comparison.In particular, in addition to being inherently dependent on the initial time t 0 , the results obtained by means of the stroboscopic Floquet approach turn out to be anisotropic, with mirrored NNN neighbors constituting complex conjugate pairs, namely t np 1,i+3 = t np * 1,i .Interestingly, we have confirmed that isotropy cannot be restored, even by considering a different initial time t 0 , akin to the findings in Ref. [24].However, as we shall explicitly compute in the following section, these 'local' features of the specific effective Hamiltonian representation are not relevant as far as the topological properties of the system are concerned.In addition, despite the different values that the NNN tunneling rates may take within each approach, it is worth noticing that they follow the predicted scaling behavior with respect to j 2 0 /ω, as it is evident from the color scale of the figure.This provides further validation of the theoretical framework.
It is worth noting that, upon numerically computing the Floquet Hamiltonian from Eq. ( 32), one can easily extract all tunneling coefficients at any order within the tight-binding expansion using Eq. ( 33).We have verified that the tunneling coefficients beyond the NNN order are indeed negligible.Therefore, in the current fast-driving scenario, the tight-binding expansion can be safely truncated at that order.Then, by using Eq.(A4), the Floquet Hamiltonian can be analytically expressed as again with h F 0 (k) ≡ 0.Here ∆ np includes also the explicit energy offset δ between sites A and B.

IV. TOPOLOGICAL PHASE DIAGRAM
We will now turn to analyze the topological properties of the system, considering both the case analyzed so far, with δ = 0, and the one in which IS is explicitly broken, δ = 0.In particular, we will draw the topological phase diagram as a function of the driving amplitudes A 1 and A 2 (for δ = 0) or as a function of A 1 and δ, by keeping A 2 fixed (for δ = 0).
The central quantity needed to discuss the topological properties of the system is the Chern number c.For a given band n, the corresponding Chern number can be written as [24,45] where Ω n (k) is the corresponding Berry curvature with ǫ n (k), |n(k) being the eigenvalues and eigenstates of the Hamiltonian h(k), while h α (k) is a shorthand notation for ∂h(k)/∂k α (α = x, y).In the present context of a two-band problem where solely the lowest band is supposed to be occupied, only the knowledge of c ≡ c 1 is required.Therefore, in Eq. (38), one can set n = 1 and m = 2.In the numerical calculations, we find convenient to tile the reciprocal space with rhomboidal cells containing the two inequivalent Dirac points, akin to the coordinate space tiling shown in Fig. 1, see Appendix C.This choice proves particularly effective for an accurate computation of the Chern number since the main contribution to the integral comes from the proximity of the Dirac points, where the Berry curvature Ω n (k) is peaked (see also Ref. [24]).In addition, we have also verified that both quasienergies ǫ n (k) (n = 1, 2) never wrap around the edges of the quasienergy first Brillouin zone [15,26], so that the denominator of Eq. ( 38) can be computed with no ambiguity.
The topological phase diagram obtained from Eq. ( 37) is shown in the left panels of Fig. 6, for the same fastdriving regime ω = 100 j 0 considered in the previous section.It is characterized by a very rich structure, in which topologically nontrivial phases with c = ±1 alternate as a consequence of the interplay between IS and TRS dynamical breaking.We recall that the current driving protocol does not exhibit IS, so that IS gets broken even in the case in which the onsite energies of undriven model are not degenerate (δ = 0), see Fig. 6a.This would not be possible in case of a monochromatic driving, like the one considered in Refs.[12,24], where an explicit breaking of IS is necessary in order to reproduce the full Haldane topological phase diagram.
Along with the phase diagram, in the right panels of Fig. 6 we show the minimal energy gap at the Dirac points, defined as ∆ǫ . This quantity is specially relevant, considering that the boundaries between the different topological phases correspond to the vanishing of the gap at one of the two inequivalent Dirac points.In particular, the thick black lines in Fig. 6b,d correspond to the regions where ∆ǫ D = 0.It is worth noticing that, though they exhibit a finite width due to numerical precision, they nicely match the phase boundaries in Figs.6a,c.This further confirms the robustness of the present results.
It is worth emphasising that here we have opted to construct the topological phase diagram based on the physical parameters characterizing the driving (A 1 , A 2 ) and the original Hamiltonian (δ), rather than relying on the conventional approach found in the literature, which employs tight-binding parameters such as the tunneling amplitude |t 1 | and its phase φ [20,41] (see also the discussion in Ref. [45]).This choice is motivated by two key considerations.Firstly, the amplitudes (A 1 , A 2 ) are the parameters that define the driving protocol, and as such they are directly accessible and tunable in any experimental implementation, making them the natural choice from the experimental point of view.Secondly, as we have seen in the previous section, the values of the tunneling parameters are inherently model-dependent, introducing complexity and potential ambiguity into the interpretation of results.
Nevertheless, it is important to remark that the two approaches for constructing the topological phase diagram are equivalent in terms of physical content.Indeed, we have verified that the topological phases predicted by the NNN truncated effective approach, as presented in Fig. 2 of Ref. [20], can be mapped onto the phase diagram in Fig. 6a.This can be conveniently demonstrated by determining the values of ∆(A 1 , A 2 ) from Eq. ( 24), along with |t (1) 1 |(A 1 , A 2 ) and φ(A 1 , A 2 ) from Eq. ( 28), and then associating the corresponding Chern numbers shown in Ref. [20] (alternatively, the Chern number can be computed as discussed in Ref. [6]).As anticipated in the introduction, this result provides a direct confirmation of the fact that the topological properties of the driven model do not depend on the specific approach employed to compute them, either the effective Hamiltonian Ĥeff or the stroboscopic Floquet Hamiltonian ĤF [t 0 ].

V. CONCLUSIONS
By means of a direct non-perturbative numerical approach, we have computed the stroboscopic Floquet Hamiltonian and the topological phase diagram of a system of non-interacting ultracold atoms in a shaken honeycomb lattice, subjected to the optimal driving of Ref. [20].Such a driving induces dynamical breaking of inversion and time-reversal symmetries, thus providing an effective way to reproduce the topology of the Haldane model, to some extent.The numerical results have been compared with the analytical prediction of the nonstroboscopic effective approach of Ref. [20], valid in the fast-driving regime, with the objective to provide a general discussion of the different perspectives for computing the effective Floquet Hamiltonian of periodically driven systems.
We have found that the topological properties are indeed gauge invariant, as expected.However, some of the tunneling parameters are model-dependent.Specifically, we confirm that the emergent NN tunneling rates are isotropic [20], meaning that the position of the Dirac points is fixed solely by the lattice geometry and is, therefore, model-independent.On the other hand, the NNN tunneling rates turn out to be anisotropic within the stroboscopic Hamiltonian approach, in contrast to the nonstroboscopic case [20].Therefore, we argue that the topological phase diagram is most effectively represented in terms of the experimentally controllable physical parameters defining the driven Hamiltonian, rather than the model-dependent emergent tunneling parameters-such as the next-nearest-neighbor tunneling amplitude |t 1 | and its phase φ [20,41,45].
Finally, we emphasize that the Floquet Hamiltonian approach represents a very valuable complementary approach to the gauge-invariant representation provided by the effective Hamiltonian of Ref. [13].In fact, in the case of a non-interacting system, it enables a direct nonperturbative calculation using simple numerical methods.In this respect, the fact that it is intrinsically dependent on the initial time t 0 is only a minor shortcoming, as all physical quantities are gauge invariant and can be computed without ambiguities.

FIG. 1 .
FIG.1.Sketch of the honeycomb lattice.The unit cell is highlighted in light green, and the sites A and B are denoted by the filled (red) and empty (white) circles, respectively.NN sites are connected by the vectors a1, a2 and a3, with j0 being the corresponding tunneling rate.The vectors b1, b2 and b3 connect the different unit cells in the Bravais lattice.Notice that there exists a specific relation between the two sets of vectors, e.g., b1 = a2 − a3; the other relations are obtained by a cyclic permutation of the indices.

FIG. 4 .
FIG. 4. (a) Onsite energy ∆ np and (b) NN tunneling amplitude t np 0 , obtained from the non-perturbative numerical approach (see text), as a function of the amplitudes A1 and A2 of the driving.Here ω = 100 j0 (fast-driving regime) and δ = 0 (no explicit breaking of IS).

FIG. 5 .
FIG. 5. Behavior of the emergent NNN tunneling rates t np 1i obtained from the non-perturbative numerical approach (see text), plotted as a function of the amplitudes A1 and A2 of the driving.The modulus and phase are presented in the left and right columns, respectively, for i = 1 (a, b), i = 2 (c, d), and i = 3 (e, f ).Here ω = 100 j0 (fast-driving regime) and δ = 0 (no explicit breaking of IS).

FIG. 6 .
FIG. 6. Topological phase diagram in the fast-driving regime ω = 100 j0.Panels (a,c) display the Chern number c obtained from the numerical non-perturbative approach: (a) Inversion-symmetric case, δ = 0, as a function of the driving amplitudes A1 and A2; (b) Explicit breaking of IS, as a function of A1 and of the bare onsite energy gap δ,for A2a/( ω) = 1 [as indicated to the dot-dashed line in panel (a)].Panels (c,d) show the corresponding behavior for the minimal energy gap at the Dirac points, ∆ǫD (see text).