Finite size effects and Hofstadter butterfly in a bosonic Mott insulator with relativistic dispersion background

TL;DR

This paper explores how strong interactions affect the Hofstadter butterfly spectrum in bosonic Mott insulators with relativistic dispersion, revealing fractal structures in various lattice geometries under synthetic gauge fields.

Contribution

It demonstrates the impact of strong on-site interactions on the Hofstadter spectrum in relativistic bosonic lattices, including finite size and temperature effects relevant for experiments.

Findings

Fractal Hofstadter spectrum appears in bosonic Mott insulators with strong interactions.

Different lattice geometries produce distinct Hofstadter butterfly patterns.

Finite temperature and size effects are analyzed within the strong coupling expansion.

Abstract

Gauge potentials with different configurations have been recently realized in the optical lattice experiments. It is remarkable that one of the simplest gauge can generate particle energy spectrum with the self-similar structure known as a Hofstadter butterfly. We investigate theoretically the impact of strong on-site interaction on such a spectrum. In particular, it is shown that the fractal structure is encoded in the quasi-particle and hole bosonic branches. A square lattice and other structures (brick-wall and staggered magnetic flux lattice) with relativistic energy dispersions which are currently accessible in the experiments are considered. Moreover, although in brick-wall and staggered flux lattices the quasi-particle densities of states looks qualitatively similar, the corresponding Hofstadter butterfly assumes different forms. In particular, we use a superposition of two different synthetic gauge fields which appears to be a generator of non-trivial phenomena in the optical lattice systems. The analysis is carried out within the strong coupling expansion method on the finite size lattices and also at finite temperatures which are relevant for the currently made experiments.

Equations

$$H=-\sum _{⟨ij⟩}{J}_{ij}{\hat{b}}_i^{†}{\hat{b}}_j+\frac{U}{2}\sum _i{\hat{b}}_i^{†}{\hat{b}}_i^{†}{\hat{b}}_i{\hat{b}}_i-\mu \sum _i{\hat{b}}_i^{†}{\hat{b}}_i,$$

$$\mathcal{Z}=\int \mathcal{D}{b}^{\ast }\mathcal{D}b{e}^{-\left(S_0+S_1\right)/\hslash },$$

$$S_0$$

$$S_1=-\sum _{⟨ij⟩}{\int }_{{}_0}^{\hslash \beta }d\tau J_{ij}{b}_i^{\ast }(\tau )b_j(\tau ).$$

$$S^{eff}$$

$$G_0^{-1}\left(\tau -{\tau }^{\prime }\right)=-{⟨b_i^{\ast }(\tau )b_i({\tau }^{\prime })⟩}_0$$

$$S^{eff}=-\sum _{ij}\sum _m{\bar{b}}_{im}\left[J_{ij}+\hslash G_0^{-1}\left(i{\omega }_m\right){\delta }_{ij}\right]b_{jm}.$$

$$S^{eff}=-\sum _m{B}_m^{†}\left[\mathbf{J}+\hslash G_0^{-1}\left(i{\omega }_n\right)\mathbf{I}\right]B_m$$

$$S^{eff}=-\sum _m{\Phi }_m^{†}\left[{\mathbf{J}}_d+\hslash G_0^{-1}\left(i{\omega }_m\right)\mathbf{I}\right]{\Phi }_m$$

$${\epsilon }_{\lambda ,\text{min}}+\hslash G_0^{-1}\left(i{\omega }_m=0\right)\mathbf{I}=0$$

$${\mathcal{G}}^{MI}\left({\omega }_m\right)=-⟨B_m^{†}{B}_m⟩=-⟨{\Phi }_m^{†}{\Phi }_m⟩=\sum _{{}_{\lambda }}{G}^{MI}\left({\omega }_m,{\epsilon }_{\lambda }\right).$$

$$G^{MI}\left({\omega }_m,{\epsilon }_{\lambda }\right)$$

$$E^{\pm }\left({\epsilon }_{\lambda }\right)$$

$$z\left({\epsilon }_{\lambda }\right)=\frac{E^{+}\left({\epsilon }_{\lambda }\right)+\mu +U}{E^{+}\left({\epsilon }_{\lambda }\right)-E^{-}\left({\epsilon }_{\lambda }\right)}.$$

$${\rho }_{MI}^{TSA}\left(\omega \right)$$

$$\rho \left(\omega \right)$$

$${\rho }_{MI}^{HSA}\left(\omega \right)$$

$$E_{p/h}({\epsilon }_{\lambda })=A({\epsilon }_{\lambda })-B({\epsilon }_{\lambda })\sin \left(\frac{\pi }{6}\mp \frac{1}{3}\arccos C({\epsilon }_{\lambda })\right),$$

$$E_t({\epsilon }_{\lambda })=A({\epsilon }_{\lambda })+B({\epsilon }_{\lambda })\cos \left(\frac{1}{3}\arccos C({\epsilon }_{\lambda })\right),$$

$$J_{ij}\to J_{ij}{e}^{i2\pi {\int }_j^i{\mathbf{A}}_L\cdot d\mathbf{l}}$$

$$\epsilon \left(k_x,k_y,\varphi =0,p/q=0\text{or}1\right)=2\left(\cos k_x+\cos k_y\right)$$

$$\epsilon \left(k_x,k_y,\varphi =0,p/q=1/2\right)=\pm 2\sqrt{{\cos }^2{k}_x+{\cos }^2{k}_y}.$$

$$\epsilon \left(k_x,k_y,\varphi =\pi ,p/q=0\right)$$

$$=\pm 2\sqrt{{\cos }^2\left(\frac{k_x+k_y}{2}\right)+{\cos }^2\left(\frac{k_x-k_y}{2}\right)},$$

$$\frac{1}{\hslash }{G}_0\left(i{\omega }_m\right)=-\sum _{n_0=0}^{\infty }\frac{\left(n_0+1\right)\left(f_{n_0+1}-f_{n_0}\right)}{i\hslash {\omega }_m-\left(E_{n_0+1}-E_{n_0}\right)},$$

$$f_k=\frac{e^{-\beta E_k}}{{\sum }_{n_0=0}^{\infty }{e}^{-\beta E_{n_0}}},$$

$$\frac{1}{\hslash }{G}_0^{TSA}\left(i{\omega }_n\right)=$$

$$\frac{1}{\hslash }{G}_0^{HSA}\left(i{\omega }_m\right)=-\frac{f_1-f_0}{i\hslash {\omega }_m-\left(E_1-E_0\right)}$$

$$-\frac{2\left(f_2-f_1\right)}{i\hslash {\omega }_m-\left(E_2-E_1\right)}-\frac{3\left(f_3-f_2\right)}{i\hslash {\omega }_m-\left(E_3-E_2\right)}.$$

$$A\left({\epsilon }_{\lambda }\right)=\frac{1}{3}\left({\epsilon }_{\lambda }+3U-3\mu -4{\epsilon }_{\lambda }{f}_3\right),$$

$$B({\epsilon }_{\lambda })=\frac{2}{3}\sqrt{3U^2+{\epsilon }_{\lambda }\left(1-4f_3\right){}^2+3{\epsilon }_{\lambda }U\left(-1+2f_1+4f_2-2f_3\right)},$$

Full text

Keywords: ultra-cold atoms, strongly correlated bosons, optical lattice, dynamics, synthetic magnetic field, vector potential, Bose-Hubbard model

I Introduction

Ultracold atoms in an optical lattice have become very clean framework for quantum simulations of condensed matter and high energy physics Bloch et al. (2008); Lewenstein et al. (2007). In particular, optical lattice systems offer the opportunity to study Abelian or non-Abelian gauge fields. Some types of such potentials have been experimentally realized in recent years Lin et al. (2009); Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Struck et al. (2012, 2013); Kennedy et al. (2015); Sugawa et al. (2016); Stuhl et al. (2015); Mancini et al. (2015). An example of the non-trivial phenomena generated by the simplest Abelian gauge (Landau configuration) is the Hofstadter butterfly Hofstadter (1976). It corresponds to the uniform magnetic field and implies fractal energy spectrum. This spectrum has been investigated in a much broader context of different gauges and lattice geometry modifications (see, e.g. Refs Goldman et al. (2009); Bermudez et al. (2010); Osterloh et al. (2005); Jaksch and Zoller (2003); Y ılmaz et al. (2015); Asboth and Alberti (2016); Zaleski and Polak (2011); Agazzi et al. (2014); Wang and Troyer (2014)). In order to get a more general theoretical description of the atoms on an optical lattice in the systems with flux attached, the effects caused by particles interactions should be taken into account Greiner et al. (2002); Bakr et al. (2010); Sherson et al. (2010); Endres et al. (2012). Therefore if we are interested in the bosonic dynamics in which the single-particle spectrum is affected by the gauge fields Lim et al. (2008, 2010); Miyake et al. (2013); Aidelsburger et al. (2011, 2013); Kennedy et al. (2015); Sajna et al. (2014); Sajna and Polak (2014); Polak and Zaleski (2013), we should take such correlations into account Powell et al. (2011, 2010); Sajna et al. (2014); Sajna and Polak (2014); Oktel et al. (2007); Umucal ılar and Oktel (2007); Möller and Cooper (2010); Huang and Wan (2011); Lim et al. (2010, 2008); Nakano et al. (2012); Qin et al. (2016); Sachdeva and Ghosh (2012); Sachdeva et al. (2010); Greschner et al. (2016); Wu et al. (2016); Ozawa et al. (2015); Vasić et al. (2015); Natu et al. (2016); Kuno et al. (2015); Toga and Yokoyama (2015); Bilitewski and Cooper (2015); Yao and Zhang (2014). It is interesting to note that recently the first experimental realization of strongly interacting bosonic atoms in a gauge field was made on the ladder lattice in the few body limit Tai et al. (2016)

Here we focus on the strongly correlated bosonic system described by the Bose Hubbard model (BHM) Gersch and Knollman (1963); Fisher et al. (1989). We use the effective field theory to investigate the quasi-particle spectrum of the bosonic Mott insulator (MI) phase. The analysis is performed for the Landau gauge in the whole range of flux strength per plaquette . So far only some chosen values of flux strength have been investigated in the BHM within superfluid (SF) and MI spectra Powell et al. (2011, 2010); Sajna et al. (2014); Sajna and Polak (2014). In particular, we analyze the appearance of the self-similar structure in a collective behavior of strongly correlated bosons.

Knowing that optical lattice patterns give us the opportunity to modify the geometry of lattice structures, we show that such modifications have a significant impact on the fractal-like structure in the MI phase (in this work as an example we analyze square and brick-wall lattices Tarruell et al. (2012); Abanin et al. (2013)). Further, we show that the tunability of self-similar spectrum of the strongly interacting bosons can be widely extended and goes beyond the geometrical type modifications of the lattice. Namely, one can consider a lattice which is already equipped by the synthetic gauge field. To show this we use one of the simplest gauge field configuration which is Abelian and also has some free control parameter, i.e. staggered flux lattice. Moreover, we deliberately choose the brick wall lattice and staggered flux lattice for the analysis. Our aim is to show that although these two types of lattices have qualitatively similar densities of states (DOS) (with relativistic dispersion), they give very different structures of the Hofstadter butterfly spectra. These differences are especially enhanced when the free parameter of the staggered flux lattice is tuned.

It is important to stress, that our calculations include the finite size effects which are present in the current experiments Stöferle et al. (2004). Moreover, to properly establish the stability of MI phase in the parameter space of BHM, the phase diagram analysis is done. We show how finite size effects and the non-trivial lattices modify the critical line which extends the previous works on this subject Lim et al. (2010, 2008); Oktel et al. (2007); Nakano et al. (2012).

It is worth adding that very recently Hofstadter spectrum was investigated in the hard-core limit of two component BHM Qin et al. (2016). The authors of this work have focused on the BHM topological properties. In contrast, we investigate a finite interaction strength which enable us to study the SF-MI phase boundary and the spectra of MI phase Bloch et al. (2008); Lewenstein et al. (2007).Besides, we also consider the effects of thermal fluctuations on the bosonic dynamics. Especially, we discuss their impact on the relativistic part of the spectrum which can be useful for future experimental setups.

The paper is organized as follows, we first describe the model and method applied (Sec. II). Next, in Sec. III, we use this method in the analysis of the self-similarity of the quasi-particle excitations for different lattices in bosonic MI phase. At the end of Sec. III we also discuss the effects of thermal excitations in the experimentally achievable range of the model parameters. Finally in Sec. IV we give a summary of our work.

II The Method

II.1 Effective action in MI phase for finite size lattice

We consider BHM with the Hamiltonian given by

[TABLE]

where $J_{ij}$, $U$, $\mu$ are the hopping energy, on-site interaction energy and chemical potential, respectively. Besides, $\hat{b}_{i}$ and ($\hat{b}_{j}^{\dagger}$) is the annihilation and creation bosonic operators at site i (the number of lattice sites is $N$). In the coherent path integral representation, the partition function of BHM takes the form

[TABLE]

[TABLE]

[TABLE]

where $b_{i}(\tau)$ is the complex field over imaginary time and $\beta$ is the inverse of temperature $1/k_{B}T$ ($k_{B}$ is the Boltzmann constant). To obtain the effective quadratic action in the Mott insulator phase, we employ the strong coupling method from Sengupta and Dupuis (2005). This method assumes perturbative treatment of $S_{1}$ in which the double Hubbard-Stratonovich transformation together with cumulant expansion are used. Then the second order effective action in $b$, $\bar{b}$ fields reads

[TABLE]

where

[TABLE]

and $\left\langle...\right\rangle_{0}=Z_{0}^{-1}\int\mathcal{D}b^{*}\mathcal{D}b\,...\,e^{-S_{0}/\hbar}$ with $Z_{0}=\int\mathcal{D}b^{*}\mathcal{D}b\,e^{-S_{0}/\hbar}$ (see also Appendix V.1). In the Matsubara frequency representation $\omega_{m}=2\pi\beta/m$ ($m\in\mathbb{Z}$), Eq. (5) yields diagonal form

[TABLE]

which can be rewritten in the matrix representation as follows

[TABLE]

where we denote $B_{m}=\left[b_{1m},b_{2m},...,b_{Nm}\right]^{T}$. Now, it is explicitly seen, that the problem of evaluating the effective action from Eq. (8) reduces only to dealing with the free particle non-diagonal part of the action $S^{eff}$ i.e. $\mathbf{J}$, because $G_{0}^{-1}\left(i\omega_{n}\right)\mathbf{I}$ is already a diagonal matrix. This is a key point of our calculations and it will be discussed more explicitly later on.

For a hopping matrix $\mathbf{J}$ in Eq. (8), one can perform unitary transformation $B_{m}=U^{\dagger}\Phi_{m}$

[TABLE]

with $\mathbf{J}_{d}=U\mathbf{J}U^{\dagger}$ where $\mathbf{J}_{d}$ is the diagonal matrix with eigenvalues $\varepsilon_{\lambda}$ , $\lambda\in\left\{1,2,...,N\right\}$ and we denote $\Phi_{n}=\left[\phi_{1m},\phi_{2m},...,\phi_{Nm}\right]^{T}$.

II.2 Finite size phase diagram at zero temperature

The phase diagram boundary is calculated from the vanishing of the second order coefficient in the effective action from Eq. (9), i.e.

[TABLE]

where this condition is met at the static limit $i\omega_{m}=0$ and for the lowest eigenvalues of $\mathbf{J}_{d}$ (which we denote by $\varepsilon_{\lambda,\textrm{min}}$). The phase boundary in terms of Eq. 10 has a mean-field character therefore we only consider the zero temperature limit in our two dimensional considerations Hohenberg (1967). Then, it is enough to use three state approximation for the local Green function $G_{0}\left(i\omega_{m}\right)$ whose details are given in Appendix V.1.

II.3 Quasi-particle density of states at zero temperature

If we are interested in the quasi-particle DOS in the MI phase, we firstly define Green’s function in the Mott insulator regime

[TABLE]

where

[TABLE]

with

[TABLE]

[TABLE]

and $G_{0}\left(i\omega_{m}\right)$ is given in Eq. (27). Then the DOS is calculated within the standard procedure

[TABLE]

where $\eta$ is the spectrum broadening parameter. It is important to stress that the form of Eqs. (12-14) corresponds to the standard form known in the Random Phase Approximation (RPA) method for the MI Green function in which eigenenergies of $\mathbf{J}_{d}$ in Eq. (9) are numerated by the wave vector $\mathbf{k}$ Sengupta and Dupuis (2005); Freericks et al. (2009); Sajna et al. (2015).

Moreover, from the equation

[TABLE]

we calculate the free single-particle DOS. Czajka et al. (2006)

II.4 Quasi-particle density of states at finite temperatures

To describe quasi-particle DOS at finite temperatures, we use the higher state approximation (HSA) for the local Green function $G_{0}(i\omega_{m})=G_{0}^{HSA}(i\omega_{m})$ (see Appendix V.1 and Eq. (28)) Sajna (2016). This form of $G_{0}$ takes into account thermal fluctuations which are observed in the periodically modulated lattice experiments Sajna (2016); Stöferle et al. (2004).

Applying similar a procedure as in Sec. II.3 but with the local Green function $G_{0}^{HSA}(i\omega_{m})$, we get

[TABLE]

where

[TABLE]

[TABLE]

and the definitions of $A(\varepsilon_{\lambda})$, $B(\varepsilon_{\lambda})$, $C(\varepsilon_{\lambda})$, $L\left(x\right)$ are given in Appendix V.2. The above calculations extend the results obtained previously for the cubic lattice with periodic boundary conditions from Ref. Sajna (2016).

In the following, we exploit the above method to study the finite-size lattices currently accesible in the optical lattice experiments, i.e. we focus on the square and brick-wall geometry and the lattice with a staggered magnetic field Tarruell et al. (2012); Abanin et al. (2013); Lim et al. (2010, 2008).

III Results

III.1 Phase diagram at zero temperature limit

Before we follow the considerations which involve the analysis of quasi-particle spectra in the MI phase, firstly we have to assess the range of Hamiltonian parameters for which the MI phase is a ground state of BHM. This will be done by investigating the phase boundary between MI and SF phase at zero temperatures by using Eq. (10). As we mentioned before we are interested in the different lattice/gauge geometries, like square and brick-wall type, and the lattice with staggered magnetic flux. Moreover, the particle hopping on these lattices will be affected by the Landau gauge which result in the self-similar structure encoded in the particle excitations. It can be achieved by proper incorporation of the hopping matrix $\mathbf{J}$ into the effective action $S^{eff}$ from Eq. (8) and by introduction of Peierls substitution to it Jiménez-García et al. (2012); Sajna et al. (2014); Kopeć and Polak (2002), i.e.

[TABLE]

where the Landau gauge is defined as $\mathbf{A}_{L}=B(0,x,0)$ (with flux per plaquette $p/q$ defined as $p/q=\eta Ba^{2}$ where $a$ is the lattice spacing and $\eta=1$ ($\eta=2$) for the square (brick-wall) lattice). Moreover, the staggered flux lattice represented by the gauge \mbox{\mathbf{A}}_{S} Lim et al. (2008, 2010), together with Landau gauge $\mathbf{A}_{L}$, are introduced to the system in superposition, i.e. \mathbf{A}_{L}+\mbox{\mathbf{A}}_{S}. Formally, the hopping around the plaquette with vertices 1, 2, 3, 4 for \mbox{\mathbf{A}}_{S} is defined as $-Je^{i\phi/4}\left(a_{2}^{\dagger}a_{1}+a_{3}^{\dagger}a_{2}+a_{4}^{\dagger}a_{3}+a_{1}^{\dagger}a_{4}\right)$ which yields $\phi$ flux, see Fig. 1 c. Moreover, for clarity, we plot the relevant lattices and gauges in Figs. 1 a and b.

The phase diagrams for the finite size lattices (30x30 sites) without uniform synthetic magnetic field ($p/q=0$) are presented in Fig. 2. On the mean-field level, the phase boundaries for the square and staggered flux lattices, have been earlier analyzed in Refs. van Oosten et al. (2001); Sajna et al. (2015); Lim et al. (2010) (phase boundary for the brick-wall lattice is the same as in the honeycomb case Li et al. (2015); Lühmann (2013)). As we see, the brick-wall and staggered flux lattices favor the MI phase in which the atoms are likely to be localized at the lattice sites. This effect could be simply accounted for the fact that both types of lattices break translational symmetry and the elementary unit cell contains two lattice sites.

Next, we focus on a more interesting situation, i.e. when the uniform synthetic magnetic field is applied in the whole relevant range of $p/q$. Moreover in our analysis, we also consider the finite size effects the lattice pattern which will be discussed in the following.

Fig. 3 gives a plot of the phase boundary at the tip of the first lobe versus flux per plaquette. Firstly, in the absence of uniform magnetic field, we compare the phase boundaries for the standard square lattice ($\phi=0$) (a), staggered flux lattice ($\phi=\pi$) (b), and brick-wall lattice (d). For the three lattices, as expected, the critical lines show different behavior. Especially we see that the critical lines for the square (Fig. 3 a) and brick-wall (Fig. 3 d) lattices are shifted with respect to each other by $p/q=1/2$ value (in the next paragraph we provide more detailed discussion of this). Moreover, it is also interesting to notice that the finite size effects are quite small for the lattices of 30x30 size for which the $p/q$ phase boundary dependence is less blurred out (compare this with 6x6 or 10x10 lattice sizes in Fig. 3 and with periodic boundary condition calculations in Ref. Oktel et al. (2007) in which standard square lattice was considered). Therefore, our calculations indicate that non-monotonous critical line behavior is quite well resolved on the 30x30 lattice size which is useful information for the future experimental setups (30x30 lattice size systems are currently realizable, e.g. see Bakr et al. (2011, 2010)).

As we mentioned above, we provide here a more detailed discussion of the critical line versus $p/q$ dependence for different staggered flux amplitudes which are parametrized by $\phi$ (see Figs. 3 a-c). In Figs. 3 a and c, the two limits of $\phi$, i.e. $\phi=0$ and $\phi=\pi$ correspond to the standard square lattice and the staggered flux lattice with maximal value of flux per plaquette, respectively. As we pointed out above both critical lines are shifted to each other by $p/q=1/2$. It can be argued on the background of the tight binding dispersions whose analytical forms for the lattice with periodic boundary conditions are known. Namely, for the case $\phi=0$ (Fig. 3 a), uniform magnetic field linearly increase from $p/q=0$ up to maximal value at $p/q=1/2$ and linearly decrease from $p/q=1/2$ up to $p/q=1$. This results in the symmetrical form of the phase boundary around the $p/q=1/2$ point. In particular at $\phi=0$ for $p/q=0$ and $p/q=1$ lattice dispersion is

[TABLE]

and at point $p/q=1/2$ the dispersion shows two branches Sajna et al. (2014)

[TABLE]

where $k_{x}$ and $k_{y}$ are wave vectors. In contrast, for $\phi=\pi$ (Fig. 3 b) the sweeping range of $p/q$ parameter is the same but at $p/q=0$ and $p/q=1$ lattice dispersion is Lim et al. (2008, 2010)

[TABLE]

which has the same lowest energy behavior as for $\phi=0$ at $p/q=1/2$ point (see Eq. (22)). Moreover, the point $p/q=1/2$ for $\phi=\pi$ and its lowest energy behavior is the same as in $\phi=0$ for $p/q=0$ and $p/q=1$ (Eq. 21). Therefore, we see that these points, and also the rest of the critical lines in Fig. 3 a) and b) correspond to each other and are shifted by $p/q=1/2$.

However the correspondence between the phase boundary shapes obtained for the standard square lattice ($\phi=0$) and the staggered flux lattice ($\phi=\pi$) for different $p/q$ parameters disappears when the phase boundary is analyzed for the intermediate values of $\phi$, i.e. for $\phi=\pi/4$, $\phi=\pi/2$, $\phi=3\pi/4$ (see, Fig. 3 c). Then the critical line changes nonlinearly between the values $\phi=0$ and $\phi=\pi$ and we cannot provide any intuitive interpretation for this situation.

Moreover, we numerically discover that at the points $p/q=1/4$ and $p/q=3/4$ all critical lines for different $\phi$ almost intersect each other (see arrows in Fig. 3 c). This intersection area becomes smaller with increasing lattice size. It is confirmed by the results presented in Fig. 4 which show that with increasing lattice size (here up to 200x200 lattice sites) the dependence of critical line at $p/q=1/4$ for different $\phi$ value is almost constant (see the scale on the vertical axis in Fig. 4). This simply shows that the tip of the first lobe in the phase diagram is intact when staggered flux of different amplitudes are applied to the system and this only happens when the system is subjected to a uniform magnetic field with $p/q=1/4$ or $p/q=3/4$ strength.

Summarizing this subsection, we have shown that a combination of the Hofstadter butterfly spectrum with a staggered flux background gives nontrivial MI-SF critical behavior. In particular, this behavior is highly nonlinear and exhibits robustness in some range of Hamiltonian’s parameters. For example, the latter effect could be of interest in quality testing of gauge field within the critical region in the experimental protocol. In further discussion we focus on the bosonic MI phase and its spectral properties, which is the main subject of concern in our study.

III.2 Hofstadter butterfly in the MI phase for different lattice geometry

Here we analyze the dependence of the MI spectrum in the whole range of the uniform synthetic magnetic field strength within the square and brick-wall lattices. The brick-wall lattice is especially interesting because the relativistic dispersion appears in its single-particle spectrum, which will be important in further discussion. For the two lattices chosen, we set the number of sites to be $30\times 30$ and as could be seen in Fig. 5, the obtained pattern for the free bosonic case ($U=0$) resembles the previous works Hofstadter (1976); Y ılmaz et al. (2015); Umucal ılar and Oktel (2007) (the calculations were made by using Eq. (16)). To realize one of the main aim of the study, we used these spectra to evaluate the strongly correlated density of states in the MI phase with unit average density ($n_{0}=1$). It was achieved by using Eq. (15) and the density plots of these spectra in terms of fluxes per plaquette ($p/q$) versus frequencies $\omega/U$ are depicted in Fig. 6.

As follows from this figure(Fig. 6) the free particle picture is changed when the interactions between bosons are taken into account $U/J\gg 1$. It is manifested by opening a characteristic MI energy gap around $\omega/U=0$. The gap divides the quasi-particle spectrum into two branches of the hole and particle like character. Moreover, we see that each of these branches also show self-similar behavior in which the Hofstadter butterfly structure is restored. In particular in the square lattice (Figs. 6 a, b, c, d) the hole (for $\omega/U<0$) and particle (for $\omega/U>0$) branches acquire the self-similar structure which corresponds to the free particle case from Fig. 6 a (the adequate correspondence in the brick-wall lattice is between the spectra in Figs. 6 e, f, g, h and Fig. 6 b). However, deep in the Mott phase in which hopping kinetic energy is weak ($J\ll U$), the Hofstadter butterfly like behavior is not well resolved due to a finite size effect and spectral broadening of the parameter $\eta$ (see Figs. 6 a, b, e, f). Therefore, future experimental setups which investigate such a strongly correlated system could partially overcome this problem by adjusting experimental parameters close to the SF-MI phase boundary, as in Figs. 6 d and h, for which quasi-particle bandwidth broadens and $J/U$ is close to the critical value (see, Fig. 4 a and d).

Finally, we can conclude that the lattice geometry and strong interactions have a significant impact on the self-similar structure of the quasi-particle MI spectra. Further, we show that the tunability of the self-similar spectrum of the strongly interacting bosons can be widely enhanced going beyond the geometrical modification of the lattice, i.e. by using staggered flux lattice.

III.3 Hofstadter butterfly in MI phase with staggered magnetic field background

In this subsection we focus on the staggered flux lattice with checkerboard symmetry Lim et al. (2008, 2010). We use Eqs. (15-16) to obtain DOS for whole range of the uniform magnetic field amplitude (see Fig. 7). The non-interacting and interacting self-similar pattern are plotted in Figs. 7 a-d and e-h, respectively. In this figures, we show the Hosftadter butterflies for different flux per plaquette $\phi$ with chosen $J/U$. In particular, we see that different kind of fractal like pattern emerges when amplitude of $\phi$ varies. This behavior is nonlinear (similar conclusion for phase boundary analysis was made in Sec. III.1). However, one can see that self similar structure for $\phi=0$ (Fig. 5 a and 6 a-d) and $\phi=\pi$ (Fig. 7 d and h) cases are shifted each other by the $p/q=1/2$ which agrees with our earlier observation in Sec. III.1.

Moreover, this is important to notice here that although the quasi-particle spectra in the MI phase for staggered flux lattice and brick-wall lattice show qualitatively similar DOS (see Fig. 8 b) and Fig. 9 a and b), they give completely different Hofstadter butterfly patterns of DOS when a uniform magnetic field is turned on. This simple argument shows that mechanism of self-similarity behavior is much more complex and DOS picture is not an adequate tool to build up an intuition about this peculiar phenomenon. Additionally, in Figs. 8 and 9 we compare some chosen quasiparticle DOS of the finite size lattice systems with exact DOS obtained for the periodic boundary conditions (see also Appendix V.3) with $p/q=0$. This shows that 30x30 lattice size quite well reproduces the shape of DOS also in the region around Dirac points i.e. at which DOS value is highly suppressed, Fig. 8 b and Fig. 9 a and b).

As follows from the above results, the tunability of the lattice in non-geometrical way, e.g. through a gauge field, can be also interesting in studying of nontrivial phenomena. Such a staggered flux lattice is a good example although it has a relatively simple Abelian structure.

To better embed the above results in the context of real experimental system, we focus further on the problem of thermal fluctuations. Namely, such fluctuations are difficult to control in the experimental protocol, especially in the strongly interacting limit Trotzky et al. (2010); Zhou and Ho (2011). Therefore, in the following, we provide an analysis of these effects within HSA approach.

III.4 Quasiparticle excitations at finite temperatures

Optical lattice experiments are not conducted at the strictly zero temperature regime Trotzky et al. (2010); Bakr et al. (2010); Sherson et al. (2010), therefore it is important to consider the effects of thermal fluctuations. I is particularly interesting to investigate a region around Dirac points in Fig. 8 b and Fig. 9 a and b. Here we focus on the whole relevant range of temperatures in the deep MI phase, i.e. we investigate temperatures up to $T/U=0.2$ which is melting point of MI phase Gerbier (2007); Polak and Kopeć (2009). It should be added that the MI phase at a finite temperature does not exist, however it is justified to discuss the MI phase below $T/U=0.2$ because MI properties can be observed up to this temperature Gerbier (2007); Polak and Kopeć (2009).

To properly catch the finite temperature regime in the MI phase, higher order energy states should be taken into account. because these states get occupied due to thermal fluctuations Sajna (2016). In the context of the lattices with relativistic dispersion, i.e. brick-wall and staggered flux lattices, we investigate this phenomena with the HSA and calculate DOS by using Eq. (17). The results are depicted in a-c for brick-wall and d-f for staggered flux lattice. In particular, the first two peaks in Figs. 10 a and d correspond to the holon and doublon excitations (from the left). The third peak in Figs. 10 a and d is a fingerprint of the triplon defects over the MI ground state. From these diagrams, two main features follow: 1) depletion of the holon and doublon excitations at the expense of triplons ones at high temperatures (see also Figs. 10 b, c and e, f), 2) the vicinity of the Dirac points are highly robust against the thermal fluctuations (see, Figs. 10 b and e). The latter observation gives the important information for the future experimental setups which shows that the interesting Dirac like physics could still be accessed even at relatively high temperatures ($T/U<0.2$).

IV Summary

We have applied the strong coupling expansion method to the BHM in the finite size lattices, which allowed us to study the energy spectrum in the arbitrary gauge fields and at finite temperatures. As an example we focused on the lattices with relativistic dispersions i.e. on the brick-wall and staggered flux ones. In particular, we have shown that strong on-site interaction of bosons does not destroy self-similar like structures in the quasi-particle spectrum, known as a Hofstadter butterfly but modifies them significantly. We have also noticed that

• quasi-particle density of states for both types of lattices studied are qualitatively similar, however they give completely different self-similar pattern, when staggered flux amplitude is tuned. This analysis was performed for a gapped phase of the Bose Hubbard model (i.e. in MI phase). Additionally, we have presented that such a remarkable fractal patterns can be only efficiently studied in the vicinity of the phase boundary, because of widening of the quasi-particle and hole energy bands for which magnetic flux dependence on energy scales is better resolved. Moreover, we have numerically shown that the phase boundary is intact over all range of staggered flux amplitudes within the uniform magnetic field at $p/q=1/4$ and $p/q=3/4$. It indicates that simple superposition of two different synthetic magnetic fields can be a generator of non-trivial phenomena in the optical lattice systems.

In this work, we have also focused on the quasi-particle excitations at finite temperatures and investigated how they are modified by thermal defects. This investigation is especially important for the experimental realization of the gauge fields in which e.g. Dirac like physics emerges (see Ref. Goldman et al. (2014) and literature therein). In particular, we have shown that the vicinities of the Dirac points in DOS are highly robust against thermal fluctuations and can be efficiently studied in the experimental setups even at relatively high temperatures (up to MI melting point $T/U=0.2$ Gerbier (2007)).

Moreover, by including finite size effects, we have simulated the lattice sizes which are currently accesible in experimental protocols Bakr et al. (2010, 2011) and we have shown that this sizes are sufficient to observe orbital magnetic field phenomena in the BHM.

Acknowledgements.

We are grateful to P. Rożek for valuable discussions. This work was supported by the National Science Centre, Poland, project no. 2014/15/N/ST2/03459 (A. S. S.).

V Appendix

V.1 Local Green function

General form of the local Green function from Eq. (6) could be rewritten in the Matsubara frequencies as a

[TABLE]

[TABLE]

where $E_{n_{0}}=-\mu n_{0}+Un_{0}\left(n_{0}-1\right)/2$ and $n_{0}$ is an integer value denoting local occupation number of bosons per site.

In this paper, we consider three (TSA) Menotti and Trivedi (2008); Ohashi et al. (2006); Freericks et al. (2009); Konabe et al. (2006); Sengupta and Dupuis (2005); van Oosten et al. (2001); Dickerscheid et al. (2003); Sajna et al. (2015, 2014) and higher state (HSA) Sajna (2016) approximations which correspond to the truncation of the sum in Eq. (25) to the three (TSA) or more (HSA) indices around some chosen value of $n_{0}$, respectively (it turns out that this is sufficient approximation to describe MI state at low temperatures Pippan et al. (2009); Ohliger and Pelster (2013)). However, at finite temperatures both of approximations gives the same zero temperature limit of $G_{0}\left(i\omega_{m}\right)$, i.e.

[TABLE]

where $f_{n_{0}}=1$ and $f_{n_{0}-1}=f_{n_{0}+1}=0$. The form of $G_{0}^{TSA}\left(i\omega_{n}\right)$ is known as TSA because it take into account three possible bosonic occupation numbers denoted by $n_{0}-1$, $n_{0}$, $n_{0}+1$. At finite temperature the situation is more involved. It turns out that higher order approximations are needed to correctly describe thermal fluctuations effects observed e.g. in periodically driven optical lattice systems Sajna (2016); Stöferle et al. (2004). Then, when the system is e.g. in MI phase with unit density ($n_{0}=1$), it is enough to choose $n_{0}=0,\,1,\,2,\,3$ which result in the following form of the local Green function within HSA method

[TABLE]

V.2 HSA coefficients of density of states

Coefficients in Eqs. (17-19) have the following forms

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $f_{m}$ are defined in Eq. (26).

V.3 Exact density of states

For depicting of the exact quasi-particle DOS in Fig. 8 and 9 we used the following lattice DOS for free bosons (we set units $J=1$):

• DOS for square lattice with dispersion $2\left(\cos k_{x}+\cos k_{y}\right)$ takes the form

[TABLE]

• DOS for brick-walll lattice with dispersion $2\sqrt{\cos k_{x}\cos k_{y}+\cos^{2}k_{x}+\frac{1}{4}}$ Abanin et al. (2013) takes the form

[TABLE]

• DOS for staggered flux lattice with dispersion $2\sqrt{2\cos\left(\frac{\phi}{2}\right)\cos k_{+}\cos k_{-}+\cos^{2}k_{+}+\cos^{2}k_{-}}$ Lim et al. (2008, 2010) where $k_{+}=\left(k_{x}+k_{y}\right)/2$, $k_{-}=\left(k_{x}-k_{y}\right)/2$ takes the

[TABLE]