Quantum phases for bosons in a magnetic lattice with a harmonic trap

The phase diagram of ultracold bosons in a 2D optical lattice subjected to a uniform magnetic flux is studied by the inhomogeneous Gutzwiller variational method. With a harmonic external potential, the density profile exhibits wedding-cake shape and the Mott phase forms a series of circular platforms with integer fillings. There are annular rings of superfluid between these platforms. We identify a lot of phases, including the unicirculation phase, the Meissner phase, the vortex phase, the reversed current phase and the vortex superfluid phase, by examining the super-currents. The circular chiral currents of the reversed current phase are opposite to other phases in annular superfluid and we analyzed the transitions between different phases.


Introduction
In recent years, the study of ultracold gases in optical lattices subjected to synthetic gauge fields has become an focus of experimental and theoretical researches [1,2]. The optical lattices are highly tunable and allow one to explore a wide range of physical parameters [3]. On the other hand, quantum gas experiments has realized artificial gauge fields, both in harmonic traps and in optical lattices [1,2,[4][5][6][7][8][9]. By using Raman lasers and a potential gradient in atomic gases, one can produce complex hopping in the bonds of optical lattices, which generate a larger effective magnetic flux acting on the neutral atoms [41,42]. This has motivated theoretical researches on the physics of strongly correlated particles in the presence of either abelian or non-abelian gauge fields which lead to various issues such as the (fractional) quantum (spin) Hall effect with bosons [10][11][12][13][14][15][16][17][18], the Hofstadter butterfly [19], the emergence of topologically protected phases [20][21][22], etc.
The synthetic gauge fields are firstly used to study the model in one dimensional optical lattice or the N-leg ladder systems [17,21,[23][24][25][26]43]. Piraud et al established the phase diagram of the strongly interacting Bose-Hubbard model defined on a two-leg ladder geometry in the presence of a homogeneous flux [23]. They demonstrated the existence of both gapless and gapped Meissner and vortex phases. Kolley et al presented the phase diagram of strongly-interacting bosons on a three-leg ladder subjected to a homogeneous magnetic flux [24]. Hügel and Paredes found that chiral ladders to be promising link systems between one-dimensional spin orbit topological (time-reversal) insulators and two-dimensional quantum Hall insulators [17]. Chiral edge currents on two-and three-leg ladders with a homogeneous flux have been observed using experimental approaches for bosons [27,28]. Atala et al reported on the observation of chiral Meissner currents in bosonic ladders [27]. Zheng et al suggested a dynamical method to observe topology of the band by Bloch oscillations [29]. They also proposed to realize a Josephson effect in momentum space by making using of the magnetic optical ladder [30]. Stuhl et al engineered an effective magnetic field in a two-dimensional lattice which brought ultracold atomic gases into the quantum Hall regime [28].
In this paper, we focus on the ground state phase diagram of bosonic atoms in a two-dimensional (2D) square lattice subjected to an artificial magnetic flux with a harmonic trap. We study the quantum phases induced by the magnetic flux by employing the inhomogeneous dynamical Gutzwiller mean-field theory (DGMFT) [31,32]. It shows that the harmonic potential force the formation of wedding-cake shape Mott platforms, between which is an annular superfluid. The vortices appear in the SF due to the magnetic flux. We find several phases in the annular superfluid region, including the unicirculation phase (U), the Meissner phase (M), the vortex phase (V) and the reversed current phase (RC), which can be identified by the chiral current. The U phase has a clockwise annular current, while the M, V, RC phase have counter-flow currents. In addition, there is a vortex superfluid phase (VSF) whose Mott platform is replaced by the SF with vortices.
The paper is organized as follows: in section 2 we introduce the magnetic Bose-Hubbard model (mBHMV) by attaching a phase to the hopping coefficient and briefly describe the inhomogeneous DGMFT in a harmonic trap. Then we investigate the influence of harmonic trap and compare the results with the model of uniform system. In section 3 we plot the phase diagram, explore novel features of the quantum phases and analyze the three different transitions. The main results are collected in section 4.

Method
We consider the Bose-Hubbard model subjected to a magnetic flux in a 2D square lattice and trapped in harmonic potential well, as shown in figure 1. The phase factor in the hopping coefficient induces an artificial magnetic flux in each cell so we call the model as the magnetic Bose-Hubbard model (mBHMV). By taking the symmetry gauge, it reads whereb j,k (b † j,k ) is the bosonic annihilation (creation) operator at 2D site { j, k} andn j,k =b † j,kb j,k is the particle number operator. The parameters J is the hopping amplitude between the nearest neighbor sites, U is the repulsive on-site interaction, and μ is the chemical potential. In this paper we take the on-site interaction U = 1 as the units. The magnetic flux per cell Φ is related to the hopping phase φ by Φ = 2φ. Here, V denotes the strength of the trap, while x j,k and y j,k are the coordinates of the sites with respect to the center of the trap.
The ground-state phase diagram without the harmonic trap (V = 0) has been studied extensively by various method [31,33,34]. In order to find the ground state of the Hamiltonian (1), we use the inhomogeneous DGMFT [31,32,35] which provides a factorized many-body variational wavefunction for the Hubbard-type Hamiltonian of the form where |n i is the n-particle Fock state at site-i (i ∈ { j, k}) and f (i) n represents the variational probability amplitude of the number state |n . n max is a cut-off of maximum site-occupation which we set n max = 7. The probability amplitudes are normalized to unity by n max n=0 |f (i) n | 2 = 1. By introducing a site-dependent complex order parameter we obtain a variational functional of the single-site Hamiltonian. The particle number at site-i is n( r i ) = n i and the average fillingn = i n( r i )/N with N the total number of lattice sites. Substituting the Schrodinger equation by the Gutzwiller wavefunction (2), one obtains the real time evolution equations [35] i d dt where ϕ i = j∈{i} ϕ j e −ikφ , with j the nearest neighbors of site i, and k depending on the locations of i, j. These equations should be solved self-consistently with the condition (3) for complex ϕ( r i ). In order to obtain the ground state of the mBHMV, we use the standard imaginary time evolution technique site-by-site. We mention that the Gutzwiller wavefunction should to be renormalized after each step since the imaginary time evolution is not unitary. The calculations are carried out for various initial configurations to reach the ground state, rather than the local minimum. At zero magnetic field without the harmonic trap, our results for the SF-MI boundary are close to that obtained by the classical mean-field method. It is not guaranteed that the mean-field description of the system could be valid for the whole parameter regimes. Nevertheless, one can expect it grasps the basic features of the mBHMV.
instead of μ i in this paper. By this way the harmonic trap V can be regarded as an effective chemical potential μ i associated to the location of the site i. The inhomogeneous chemical potential breaks the periodic and translational symmetry of the optical lattice and has significant impact on mBHMV.
There are many studies on the boson lattice systems in a magnetic field [33,34,36]. The superfluid regions in Bose-Hubbard model are divided into particle-superfluid (p-SF) and hole-superfluid (h-SF), which have been introduced by our former work [33], as shown in figure 2. The particle density of p-SF is slightly smaller in the vortices core, while the particle density of h-SF is slightly larger in the vortices core. On the other hand, for Bose-Hubbard model in the presence of a harmonic potential without magnetic flux, superfluid and Mott-insulating phases coexist and a wedding cake structure [37][38][39] is formed in the particle density n( r ) picture. It is similar for the magnetic Bose-Hubbard model in a harmonic trap. We show the order parameter |ϕ( r )| and the particle density n( r ) of the mBHMV with Φ = π/2 4 , which contains a special 'wedding cake' in figure 3. The distribution of n( r ) can be regarded as 'wedding cake with candles'. The special superfluid regions filled with vortices ('candles') are surrounded by annular Mott regions. The phase crossed by the red arrow in figure 2, including p-SF, h-SF and Mott insulator (MI), can be found with proper parameters. To identify these phases, we plotted some cross sections of the particle density as figure 3(c), where we can find the p-SF phase (the red solid line with  some hollows, at y = 29 in figure 3(b)), the h-phase (the yellow dotted-dash line with some bulges at y = 11, the green line at y = 16) and the Mott insulator phase (the flat parts of all lines). The hollows (bulges) are particle-(hole-) type vortices.

Phase diagram
We present the ground-state phase diagram of BHM in a harmonic trap with magnetic field as a function of hopping coefficient J and magnetic flux Φ in figure 4. The phase diagram is displayed for Φ ∈ [0, π] while the whole period [0, 2π] are axisymmetry by the equation (1). Here we set equivalent chemical potential μ = 1.5(center) → 0.5(edge). The number and the width of SF ring is influenced by the equivalent chemical potential μ, the magnetic flux Φ and the hopping coefficient J. The ends of the red arrow in figure 2 are affected by μ, and the X-coordinate position of the arrow is affected by J. The different phases are primarily characterized by the super-current configurations in the superfluid ring while the MI rings are trivial. The super-current density j ≡ (j x , j y ) is calculated by To measure the super-current of different phases, we introduce the polar-current which presents clearly magnitude and direction of the super-current, where j is the super-current density, r is the radial unit vector and z is the Z-axis unit vector. When the super-current at ith ring flow anticlockwise, j ρ (i) > 0, otherwise j ρ (i) < 0. According to the different types of super-current, we have catalogued many quantum phases: (i) the unicirculation phase (U), which shows a toroidal current; (ii) the Meissner phase (M), which contains two

Unicirculation phase
The super-current (including the magnified part), which contain an only clockwise toroidal current, the polar current of the unicirculation phase, the order parameter |ϕ( r )| and the particle density n( r ) are shown in figure 5. The particle density contains a double-layer Mott platform and a superfluid ring. And the clockwise super current is formed on the superfluid ring. There are two factors to affect the toroidal current: the smaller hopping coefficient J and the magnetic field along the positive Z-axis where B = (0, 0, 2B). And we will analyze these factors in the section 3.6. Figure 5(b) plot the magnified part of figure 5(a) to observe the super-current clearly. The trough of polar current, shown in figure 5(c), corresponds to the direction of the clockwise current. Figure 6 shows the super-current, the polar current, the order parameter |ϕ( r )| and the particle density n( r ) of the Meissner phase. In figure 6(a), there are two opposite toroidal currents in the SF ring. The outer current is clockwise, while the inner current is counterclockwise. We provide an enlarged picture in figure 6(b). The two toroidal currents that appeared in the edges of the same superfluid ring are made up of a series of small semicircular orbits. The equivalent particles execute as their cyclotron orbits, then impact the Mott region wall and bounce back to the next orbit. It is shown as a schematic picture with the red arrows in figure 6(b). The direction of the current is affected by the orbits which can be seen as incomplete vortices. The polar-current of the Meissner phase, as shown in figure 6(c), is different from the unicirculation phase. The peak and trough of the polar-current respectively are related to the outer and the inner super-currents.

Vortex phase
As increasing the magnetic flux Φ or hopping coefficient J beyond a critical value, the system enters the vortex phase, resulting in a special super-current configuration filled with vortices in the SF ring. These vortices interact repulsively with each other, and unevenly distributed in the superfluid region [25], as   figure 6(c)). That means the total outer (inner) current is also clockwise (counterclockwise).

Reversed current phase
When the magnetic flux Φ exceeds about 0.84π, a reversal of the currents in the annular SF region is observed. The super-currents, the order parameter and the particle density of the RC phase are shown in

Vortex superfluid phase
With the increase of hopping J, the red arrow in figure 2 moves right and the superfluid region expands, leading to the disappearance of the central Mott domain. This phase, that contains a large superfluid filled with many vortices, cannot be classified by the super-current and the polar current, and we call it vortex superfluid phase (VSF) as shown in figure 9. There are more vortices in the margin of the whole superfluid region, which is a deeper superfluid region with a suitable equivalent chemical potential μ. The vortex superfluid phase does not have rigorous boundaries with other phases marked with dashed lines in figure 4. We think the boundaries between VSF and other phases are some smooth crossover areas. The polar-current is up and down because of the wide range of the circular superfluid region caused by the hopping J and the harmonic trap V.

Discussion
The phase diagram (figure 4) contains three important transitions: Meissner-vortex transition, unicirculation-Meissner transition and the chiral polar current reversal transition. For Meissner-vortex transition, the hopping coefficient J (corresponding to the width of the superfluid ring) and the magnetic flux Φ (corresponding to the size of vortex) affect the forms of the super-current in the SF regions of different phases together. For the Meissner phase, the edges of SF (near to the MI) have two reverse chiral super-currents. There will be unidirectional current on the SF side, which can be interpreted as the skipping cyclotron orbits as shown in figure 6(b) [40]. The direction of chiral current depends on the position of the skipping cyclotron orbits. When the skipping cyclotron orbits are on the outside of the SF, the super-current is anticlockwise; when the orbits are on the inside, the super-current is clockwise. The vortex phase prefers a wide superfluid ring and smaller vortices compared to the Meissner phase because of the repulsive interactions between vortices. On the other hand, a bigger magnetic flux leads to a smaller size of vortices. If there is enough space for particles, the large magnetic flux will destroy the unidirectional chiral currents and the Meissner phase, leading to a vortex phase.
Then we compare the unicirculation current phase and the Meissner phase. Based on the same reason, there are repulsive interactions between two opposite currents in M phase. So only one current can be left as  In order to show the current reversal transition, we introduce the difference of the annular current between the outer-ring and the inner-ring of the polar-currents j ρ with different hopping coefficient J = 0.015, 0.02, 0.03, 0.05 in figure 10. Hence, Δj ρ < 0 means reversed current while Δj ρ > 0 means normal currents. As the effective magnetic flux Φ increases, the reversed current appears. The current reversal is obtained mainly because of the enlarged unit cell and the particles with a large wavelength experiencing the effective magnetic flux [24,25]. That means there is an M-times enlarged block of the unit cell whose effective magnetic flux is MΦ leading to a reversed current when π/2 < Φ < π (set M = 2 here) [25]. The reversed current transition is caused by the competition of the magnetic flux between the unit cell (Φ) and the enlarged block of the unit cell (2Φ).

Summary
We have investigated the quantum phases of ultracold bosonic atoms in a 2D optical lattice subjected to a homogeneous magnetic flux with a harmonic potential by the DGMFT. The harmonic potential shapes the 'wedding-cake', while the magnetic flux adds the 'candles' (vortices). We plot the phase diagram of the mBHMV with a harmonic trap, including the unicirculation phase, the Meissner phase, the vortex phase, the reversed current phase and the vortex superfluid phase. The different phases are distinguished by the patterns of chiral currents. The U phase has a clockwise annular current, while the M, V, RC phase have two opposite currents.