Magnetic Shell Structure of 2D-Trapped Fermi Gases in the Flat-Band Lieb Lattices

We investigate the magnetic shell structure of repulsively interacting two-component Fermi gases trapped in a two-dimensional harmonic potential and loaded on the optical Lieb lattices. We employ the real-space dynamical mean-field theory (R-DMFT) to explicitly consider the trap potential in a self-consistent way. Computing the profiles of particle density and local magnetization across the lattice sites in the trap, we find that the incompressible core with ferrimagnetic ordering appears with the density plateau at the trap center, which is surrounded by the shell of the normal metallic phase. We examine the incompressibility of the core by adding more particles and creating the higher spin-population imbalance. While the core area expands from the outer shell with added particles and increased polarization, the excess particles are prohibited from going inside the core, and thus the density plateau is unchanged at the half-filling with the same magnetic ordering. In addition, we find that the feature of the phase separation differs with the sublattices, where the interstitial sites causing the flat band dispersion shows the signature of the abrupt transition in the density and magnetization at the boundary between the core and surrounding shells.


Introduction
Ultracold Fermi gases have attracted much attention because of its potential as an ideal test bed of condensed matter phenomena that can simulate or emulate a designed quantum many-body Hamiltonian in highly controllable environments [1][2][3][4][5]. The atomic Fermi gases confined in an optical dipole trap realize an artificial system at a temperature considerably lower than the Fermi temperature, being driven to the quantum degenerate regime [6][7][8]. The pseudospin components of the fermionic atom are given by its hyperfine states of which two are typically used in the popular systems of 6 Li or 40 K atoms, providing the up-spin and down-spin species of fermions. The scattering between these two-component fermionic atoms can be controlled by using the Feshbach resonance [9], providing a way to tune the particle-particle interaction to be attractive or repulsive. Loaded on optical lattices, fermionic atoms with tunable interactions can provide a realization of the Fermi Hubbard Hamiltonian [10]. In particular, the seminal experiments [11,12] realizing the repulsive Hubbard model demonstrated the Mott-insulating phase of the Fermi gases in cubic optical lattices in a trap by examining the bulk compressibility of the atomic cloud. Theoretical calculations using the real-space dynamical mean-field theory [13] found that the Mott-insulating core appeared with the plateau in the density profile in the three-dimensional trap.
The range of geometries of the optical lattices realized in ultracold quantum gas experiments has become much broader in recent years. The triangular, hexagonal, Kagome, and other non-standard types of lattices, including systems in dimensional crossover, are experimentally realizable [14][15][16][17][18][19][20][21][22][23][24][25][26][27], providing new opportunities to explore the physics of strongly correlated particles within the Hubbard Hamiltonian, realized with ultracold Fermi gases. One of the most recent realizations is the Lieb lattice system that is the two-dimensional lattices decorated with an interstitial site residing at the center of each bond connecting the two nearest neighbors of the regular square lattices. While it has been elusive in conventional solid-state materials, structures like the Lieb lattices that support flat-band dispersion have attracted steady interest because of their peculiar properties related to the flat-band dispersion, which introduces a large degree of kinetic frustrations [28]. Perhaps the most well-known example is flat-band ferromagnetism [29][30][31][32], which was rigorously proven in the presence of interactions in the Lieb lattices. In the context of ultracold gases, a possibility of condensation with kinetic frustrations in flat bands has been theoretically studied [33]. Very recently, there has been rapid progress in the experimental realization of Lieb lattices in artificial quantum systems. The experimental realizations of the optical Lieb lattice with ultracold Bose gases have been reported [34,35]. Also, the photonic, polaritonic, and electronic realizations of the Lieb lattices and the one-dimensional analogs of the Lieb lattices have been reported [36][37][38][39][40][41].
In this paper, we investigate the influences of the trap confinement of the ultracold gas system on the magnetic phase structure of repulsively interacting two-component fermions loaded on optical Lieb lattices within the description of the repulsive Hubbard model in an external harmonic potential. While the attractive and repulsive Hubbard models in the Lieb lattices and other flat-band lattice geometries have been an active subject of the flat-band superconductivity [42][43][44][45][46][47][48][49][50][51], magnetism [52][53][54][55][56][57][58][59][60][61][62][63][64], and topological states [65][66][67][68][69][70][71], the inhomogeneity due to the essential presence of the trap potentials and possible formation of spatial phase structures have not been considered yet in the context of the ultracold gases. In our numerical calculations, we explicitly consider the trap potentials by employing the real-space dynamical mean-field theory (R-DMFT) that includes the local quantum fluctuations beyond the mean-field theory and the site-dependent self-energies to present the density and polarization profiles in the trap in a self-consistent way. At the half-filling in the absence of the trap, a previous DMFT study with the numerical renormalization group [53] found a ferrimagnetic ordering where the staggered magnetization appeared while the system was ferromagnetic. Namely, while total up-spin and down-spin populations are not equal in an entire system, which is thus a ferromagnet, the majority spin species at a site alternates between up-and down-spin species across the lattice sites. A recent DMFT study with quantum Monte Carlo calculations [63] also found that at a low temperature, the double occupancy sharply drops at one of the sublattices because of the localization in the flat-band at the half-filling.
We find that the system in the trap exhibits a double-shell structure composed of the incompressible core of the ferrimagnetic ordering and the surrounding outer shell of the normal metallic phase at a low temperature. The core at the trap center has a finite radius, inside which the particle density is uniform with half-filled lattices, showing the plateau in the density profile. In the magnetization profile, the core shows the staggered pattern, while one of the spin species is partly depleted, verifying the ferrimagnetic ordering previously expected at the half-filling. We also examine the incompressibility of the core by adding particles and increasing the total polarization. While the core area expands with more particles being included in the trap, the excess particles do not increase the density of the core, so that the plateau in the density profile does not change. In addition, we find that the signature of the phase separation is different between the two sublattices where the discontinuous change of the density is evident only at the interstitial sites at the shell boundary between the areas of the ferrimagnetic, incompressible, and normal metal phases. This paper is organized as follows. In Section 2, we briefly review the real-space dynamical mean-field theory and provide the numerical details. In Section 3, we present our numerical results of the density and polarization profiles across lattices in the trap. The explicit real-space results are compared with the estimates based on the local density approximation. The incompressibility of the ferrimagnetic core in the trap is presented with a change of the particle number and polarization. In Section 4, we summarize our results and provide our conclusions and outlook.

Real-Space Dynamical Mean-Field Theory
We describe the trapped two-component Fermi gases loaded on optical Lieb lattices with repulsive onsite interactions by the Hubbard Hamiltonian in an external trap potential as: whereĉ † i,σ andĉ i,σ are the creation and annihilation operators of the fermionic particle with a pseudospin component σ at site i, respectively, andn i,σ is the corresponding number operator. The summation with i, j runs over the nearest neighbor pairs of the sites i and j. The physical constantsh and k B are set to be unity, and the temperature is in the energy unit of t, which is omitted as we set t = 1 in our notation hereafter. The harmonic trap potential is given as V(x, y) = 1 2 V 0 (x 2 + y 2 ) with isotropic trap strength V 0 . The particle number N = N ↑ + N ↓ = ∑ i,σ n i,σ and the total polarization P = (N ↑ − N ↓ )/N can be controlled by tuning the chemical potential µ σ . We employ the real-space dynamical mean-field theory to solve this Hamiltonian by explicitly considering the inhomogeneous trap potential.
The dynamical mean-field theory (DMFT) is a well-established beyond-mean-field method taking full account of local quantum fluctuations, which allows us to study the Mott metal-insulator transition in a unified framework [72]. However, there is a problem in the direct application to the ultracold Fermi gas systems, which originates from the fact that the trap potential and the consequent inhomogeneity is essential to the realization of the system. The easiest but approximate extension of the conventional theory constructed in a homogeneous system is to use the local density approximation (LDA) by assuming local homogeneity. In the LDA scheme, the observable that we want to know at position r in the trap is approximated from the solution of the translation invariant system by applying local chemical potential µ = µ − V(r), which is shifted by the trap potential at r. While the LDA is convenient since it does not require any modification of the traditional solid-state model, its applicability to the description of the spatial-phase separation is often questioned because of the large correlation length at phase transitions and a possibility of spatially inhomogeneous phases, such as the Fulde-Ferrell-Larkin-Ovchinnikov state of exotic superconductivity [73]. The real-space DMFT (R-DMFT) for trapped ultracold Fermi gases [13,74,75] extends the idea of DMFT to directly include the site-dependent self-energies, while it still keeps the local quantum fluctuations and self-consistency for all lattices in the trap.
The DMFT procedures map the original lattice problem onto the local quantum impurity problem. After solving the impurity problem, the local self-energy is given back to the lattice problem, and the procedures are repeated until the self-consistency in the local self-energy or Green's function is gained. In comparison, the R-DMFT has the original lattice problem mapped onto the local impurity problem for every single site, and thus we have the impurity problems to be solved as many as the number of lattice sites in the system. The impurity problems in the R-DMFT can be solved separately for each site, but the resulting site-dependent self-energies are sent back together into the original lattice problem to proceed for the next iteration. These procedures are formulated as follows.
Let us begin with the diagonal matrix of the local self-energies [Σ(σ; iω n )] ij = δ ij Σ j (σ; iω n ), where ω n is the fermionic Matsubara frequency, and δ ij is the Kronecker delta. The site-dependent Σ j (σ; iω n ) is an initial guess at the first iteration. Then, given the self-energy matrix Σ(σ; iω n ), the Dyson equation for the lattice Green's function G(σ; iω) can be written at each Matsubara frequency ω n as where G 0 is the lattice Green's function of the non-interacting system, calculated as for the trap potential V(x i , x j ) and the hopping part [H 0 ] ij between sites i and j in the original Hubbard Hamiltonian. Numerically solving Equation (2) to obtain G(σ; iω n ) requires an efficient linear solver, which is the first time-consuming part in the R-DMFT iteration. For the Lieb lattice systems in a trap, we consider up to 2700 sites which corresponds to a 2700 × 2700 matrix of G for each of the up-and down-spin components at a given frequency. The best parallelization efficiency in solving Equation (2) can be obtained if distributed over the Matsubara frequency domain.
Having the lattice Green's function G calculated from the lattice Dyson equation, one can generate the site-dependent Weiss mean-field from the local Dyson equation as: The Anderson impurity model is constructed at each site by using G (0) i which describes the quantum bath coupled to the impurity site. The numerical details of the description of the bath orbitals depend on the particular method to be used to solve the impurity model. There are a variety of numerical methods available for the impurity model in the DMFT, which include quantum Monte Carlo methods [76,77], numerical renormalization group [78], density matrix renormalization group [79], and exact diagonalization [80,81]. Here, we employ the exact diagonalization, as while it includes an approximation step in the local Hilbert space projection, it is fast enough to cope with the numerical difficulty in handling a large number of the impurity problems generated in the R-DMFT procedures [75].
The exact diagonalization scheme describes the bath Green's function by using a finite number of orbitals [72]. The corresponding single-impurity Anderson Hamiltonian is written as: (5) whereâ † andâ are the creation and annihilation operators of the bath orbitals, and the parameters E α,σ and V α,σ represent the energy of the bath orbitals and the coupling to the impurity site. These parameters are determined from the projection of the non-interacting impurity Green's function to G (0) i (σ; iω n ) given by the local Dyson equation. This projection is done in practice by numerically minimizing the objective function: where the factor 1/ω n is introduced to give more weight to the low-frequency part. Once the parameters are determined, the exact diagonalization of H i,n s straightforwardly provides the local self-energy Σ i (σ; iω n ) by using either the Lehmann representation at finite temperature, or the Lanczos method at zero temperature [72,80]. Constructing and solving the impurity Hamiltonian for every site is the second time-consuming part of the R-DMFT interaction. The computational time is reduced by parallelizing the calculations in the lattice site domain. The local self-energies collected after sweeping over all sites are sent back to the lattice Dyson solver in Equation (2). This closes one R-DMFT iteration that goes on until the local self-energies are converged. Figure 1a displays the structure of the Lieb lattices, where the unit cell consists of three sublattice sites marked as A, B, and C. We define L × L lattices as the L × L tiles of the unit cell, which includes N = 3L 2 sites in total. The origin (x, y) = (0, 0) of the harmonic trap potential V(x, y) is located at the central A site of the lattices. We set the spacing between the nearest neighboring sites to be one. We fix the temperature at T = 0.05 throughout our R-DMFT calculations in this paper. We use n s = 5 for the number of the bath orbitals to construct the Anderson impurity model, which we solve by using the exact diagonalization. The frequency grids are set with 1000 Matsubara frequencies.

Finite-Temperature R-DMFT in the Lieb Lattices
The temperature T = 0.05 is low enough to exhibit the magnetic ordering as expected in the previous DMFT calculations at the half-filling in the Lieb lattice without a trap [53,63]. Here, we mainly consider the interaction strength U = 2 where we do not see other orders, such as the strip order or the charge and spin density waves [59,63] which may occur with finite doping at lower temperatures or at stronger interactions. Note that here we reserve "R-DMFT" exclusively for the calculations with the trap potential, while we use the term "DMFT" for the translation invariant system with the three sublattices which brings the Brillouin zone summation into Equation (2). Figure 2b presents our DMFT calculations of S z ≡ 1 2 ∑ α=A,B,C (n α,↑ − n α,↓ ), where n α,σ = n α,σ , as a function of temperature T at the half-filling, verifying that the system at T = 0.05 exhibits S z = 1/2 in agreement with the ground-state DMFT calculation of the flat-band ferromagnetism [53]. Figure 2 displays the R-DMFT calculations of the particle density and polarization in the Lieb lattices of the 30 × 30 supercell structure confined in the trap potential of strength V 0 = 0.01. The chemical potential is fixed at µ ↑ = µ ↓ = 1, which would give the half filling if the trap potential is absent. In the presence of the trap, we find that there exists a shell structure in the density and polarization profiles where the polarized core at the trap center is surrounded by the unpolarized outer shell in the trap edge. The trap used in Figure 2 contains N ≈ 900 particles with the total polarization P = (N ↑ − N ↓ )/(N ↑ + N ↓ ) ≈ 0.07. Figure 2a shows the up(down)-spin-majority sites in the red(blue) color, and the spin-density-balanced sites in the gray. The size of the marker scales with the filling factor of the particles at the lattice site, indicating that the core shell has an almost uniform density of particles, while the density in the outer shell decays as it goes to the trap edge.

Magnetic Shell Structure and Phase Separation in the Trap
The profiles of the density n ≡ n ↑ + n ↓ and magnetization m ≡ n ↑ − n ↓ along the line of y = 0 across the trap characterize the different phases of the core and outer shells. In the density profile shown in Figure 2b, the density plateau at n = 1 defines the core region at the trap center, whose incompressibility will be further examined below. In the surrounding outer shell, the densities are quite different between the sublattices A and B/C, and the sublattice B/C also shows an abrupt change in the density at the boundary between the core and surrounding outer shells.
On the other hand, in the magnetization profile, the ferrimagnetic order appears in the core region, which is in good agreement with the previous DMFT result at the half-filling in the absence of the trap [53]. The sign of the magnetization m(x, y = 0) spatially alternates between positive (up-spin majority) and negative (down-spin majority) along the sublattices A and B/C, as shown in Figure 2c. While this staggered magnetization is similar to the antiferromagnetism, it is, in fact, a ferromagnet since there are more up-spin particles than down-spin particles in the core. This magnetization profile also shows the plateau in the core region, which is abruptly dropped to zero as it enters into the region of the outer normal metal shell of the trap edge. The phase separation between the magnetic incompressible core and the normal metallic surrounding appears as a sharp boundary between the shells. The rather abrupt change is observed in both of the density and magnetization profiles, particularly at the sublattice B/C, and there is an effect of the finite temperature and the influence of the finite-size system, which may be responsible for the thin transition area at the boundary. Other types of orders, such as the strip order and the charge and spin density waves predicted in References [59,63] may appear at the shoulder of the trap away from the half-filling, but we do not see them in the systems that we have examined at U = 2 and T = 0.05. However, we have experienced a numerical instability problem at the boundary sites at lower temperatures and stronger onsite interactions, suggesting a possibility that other orders may exist indeed in the ground state in the strong interaction regime.

Verification of the Local Density Approximation
One of the implications of the R-DMFT calculation for a system in a trap is that it can examine the validity of the local density approximation (LDA). Within the LDA, a local quantity of the system in an external potential can be treated as the one computed in the translation invariant system at the chemical potential shifted by the value of the trap potential at the particular location where the quantity is estimated. Here, in particular, we compare the radius of the ferrimagnetic core between the estimate based on the LDA and the shell boundary observed in the R-DMFT calculation that explicitly includes the trap potential.
In order to obtain the LDA estimate of the core radius, the DMFT calculations are given for the average density and total spin of the three sublattices as a function of chemical potential µ = µ ↑ = µ ↓ , which are shown in Figure 3. It indicates that the transition to the normal metallic phase occurs at the chemical potential µ c ≈ 0.62. Then, the radius of the core shell can be estimated within the LDA by solving the equation, where µ 0 = U/2 = 1 indicates the chemical potential used in the R-DMFT calculations, and V(R c ) = 1 2 V 0 R 2 c is the trap potential at the radius R c . The LDA estimate of the radius R c depends only on the trap strength V 0 , once µ c is determined. Table 1

Compressibility Tests of the Core by Adding Particles and Increasing Polarization
In the sense of LDA, the density plateau with n = 1 of the core shell indicates that ∂n ∂µ = 0, suggesting zero compressibility. Energetically, this is due to the spectral gap [53] within which the change of chemical potential does not induce the change of the density. Thus, putting a particle into the trap, it is energetically difficult to go into the core shell in the trap center. At best, it would go to the boundary of the outer metallic shell and might participate in the core from the outermost part if it is energetically favored than being in the normal metallic phase.
By using the R-DMFT calculations with the explicit inclusion of the trap potential, we examine more directly the effects of adding particles and increasing spin-population imbalance on the shell structure in the trap. Since R-DMFT is based on the formulation in the grand canonical ensemble, we control the spin-species-averaged chemical potential µ 0 and the Zeeman field h that determines the spin-dependent chemical potentials, µ ↑ = µ 0 + h and µ ↓ = µ 0 − h, to add more particles and create the higher spin-population imbalance. Figure 4 presents the effects of increasing µ 0 at h = 0 on the density and magnetization profiles along the axis across the trap center. While the core area expands with more particles being added, the density plateau of n = 1 does not change in the core shell. This indicates that the added particles do not go to the interior of the core to raise the density at the core center, suggesting that the ferrimagnetic core is not compressible. The added particles go to the shell boundary and then enlarge the ferrimagnetic core from the outermost region.
As the core area at constant density n = 1 gets larger with added particles, the spin population imbalance in the entire system also increases because the ferrimagnetic phase of the core has the fixed magnetization per site and per unit cell. Increasing µ 0 with h = 0 automatically leads to the higher polarization, accordingly, as the number of particles gets increased in the core. The next question would be, what would happen if we changed P to make the spin-population imbalance even larger than the total magnetization, that can fit in the core without breaking the ferrimagnetic phase? While it is numerically very time-consuming to tune µ 0 and h to fix the total particle number and spin polarization separately, here we simply increase the Zeeman field h at a fixed µ 0 to increase the polarization, which, in this case, accompanies the increase of the particle number as well. It turns out that the excess polarization goes to the metallic surrounding to make the outer shell polarized. Figure 5 presents the effects of increasing the Zeeman field h with the average chemical potential being fixed at µ 0 = 1. We observe that the ferrimagnetic core expands with increasing h, and the surrounding outer shell gets polarized at a finite h since it is in the paramagnetic metallic phase. In the tests with the several values of h, it turns out that the core of the ferrimagnetic phase is quite robust against the increase of the Zeeman field, implying that it is magnetically incompressible as well. Thus, adding spin-polarized particles of a particular spin species would not go inside of the core to increase the local density and magnetization at the trap center.
Finally, we discuss the feature of the phase transition appearing in this trapped system and its sublattice-type dependence. In a harmonic trap, the phase of the quantum state depends on the position if the phase transition is driven by changing the density of particles or the local chemical potential in the sense of the LDA. In the system of the Lieb lattices in a trap, while the phases of the core and surrounding outer shells are distinguished by the staggered magnetic ordering and the density plateau, the discontinuity characterizing the phase separation appears differently in the sublattice sites A and B/C. The jump in density is evident only at sites B/C at the boundary between the core and surrounding outer shell. Also in the magnetization profiles, the jump at the boundary is much clearer at the sites B/C, as shown in Figure 5. These differences between the sublattice sites A and B/C may originate from the different contribution of the flat-band dispersion, as previously discovered by the zero-temperature DMFT calculation in Reference [53] where the flat-band feature appears in the local density of states only at the sublattice B/C.

Summary and Conclusions
In this paper, we have investigated the flat-band ferromagnetism in the Hubbard model of repulsively interacting Fermi gases confined in a harmonic trap and loaded on the optical Lieb lattices by using the real-space dynamical mean-field theory. We have observed the spatial phase separation in the trap where the incompressible core with the ferrimagnetic ordering is surrounded by the outer shell of the normal metallic phase. The core is also identified by the density plateau of the half-filled lattice sites within the critical radius from the trap center, which is in good agreement with the estimate based on the local density approximation applied to the translation invariant system without a trap. We have also examined the responses of the phase structure of the core and surrounding outer shells to the changes of the particle number and the spin-population imbalance. We have found that the added particle and increased polarization enlarges the area of the core shell, though while keeping the density plateau and ferrimagnetic ordering unchanged-indicating that the core is incompressible. On the other hand, the surrounding outer shell of the metallic phase is polarized by the excess particles with the spin-population imbalance that are not contained by the core of the fixed density and magnetization. In the trap, the phase transition between the incompressible ferrimagnetic phase and the metallic phase is driven by the position-dependent local chemical potential and making a phase boundary between the core and outer shells. In the Lieb lattices, it turns out that the feature of the transition differs locally with the type of the sublattices at the boundary. The interstitial sites show a sharp jump in the density and magnetization profiles, which may originate from the very high local density of states due to the contribution of the flat-band dispersion. We have considered a relatively weak repulsive interaction in this study. An interesting subject to focus on for future study may be the investigation of a trapped Lieb-lattice system in intermediate and strong interaction regimes to see how the expected strip order or charge and spin density waves could be stabilized with trapped ultracold Fermi gases in the decorated optical lattices.