{'fraction_non_alphanumeric': 0.04728305175273381, 'fraction_numerical': 0.01146937283895089, 'mean_word_length': 4.3785908678560626, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Filling-driven Mott transition in SU(N) Hubbard models

We study the filling-driven Mott transition involving the metallic and paramagnetic insulating phases in SU(N) Fermi-Hubbard models, using dynamical mean-field theory (DMFT) and the numerical renormalization group (NRG) as impurity solver. The compressibility shows a striking temperature dependence: near the critical temperature, it is strongly enhanced in the metallic phase close to the insulating phase. We demonstrate that this compressibility enhancement is associated with the thermal suppression of the quasiparticle peak in the local spectral functions. We also explain that the asymmetric shape of the quasiparticle peak originates from the asymmetry in the underlying doublon-holon dynamics.


I. INTRODUCTION
The Mott insulator transition [1], as an ubiquitous phenomenon in strongly correlated systems, appears at the boundary of various phases. One important example is the transition between metal and paramagnetic insulator in the Fermi-Hubbard model, where no symmetry (e.g., spin, discrete translational invariance) is broken across the transition. Indeed, this paramagnetic Mott transition originates purely from the competition between Coulomb interaction and kinetic energy where Mott insulators allow only integer number of particles at each site, without charge fluctuations, away from completely empty or full occupation.
Consequently, the paramagnetic Mott transition can be induced by changing the chemical potential, i.e., the filling [2][3][4][5][6][7]. This filling-driven Mott transition has been studied in realizations of the Hubbard models using ultracold atoms [8][9][10][11][12], in which a harmonic confinement potential is applied to the optical lattice to impose a lattice boundary. As the confinement potential is not uniform, the filling of each site varies from site to site and so different phases can appear in different regions within a single trap of atoms.
However, the filling-driven paramagnetic Mott transition has not been observed yet in ultracold atom experiments, since this requires sufficiently low temperatures relative to the Fermi energy. While the paramagnetic Mott insulator has been observed at temperature of about 20% of the Fermi energy [8][9][10][11][12], its evolution towards a metal occurs not via a transition but crossover at such elevated temperature. An actual transition occurs only below the critical temperature which is only a few percent of the Fermi energy [13][14][15].
Recently, it was demonstrated that the one-band Hubbard model with flat potential profile can be cooled down to host the anti-ferromagnetic phase [16]. Thus one may expect that the Hubbard model with non-uniform potential also can be brought below the critical temperature of the paramagnetic Mott transition. If so, how can one discriminate between the transition and the crossover by using quantities accessible in experiments where temperature is not even directly measurable?
In this work, we study compressibility as a function of the occupation number, which has been measured in ultracold atom experiments [11,12], in the SU(N ) Hubbard model (2 ≤ N ≤ 5) with strong Coulomb interaction. We demonstrate that it exhibits distinct behaviors depending on whether the temperature T is below, above, or near the critical temperature T c . To study the paramagnetic phases of the multi-flavor model for arbitrary interaction strength, filling, and temperature, we use dynamical mean-field theory (DMFT) [17,18] and the numerical renormalization group (NRG) [19,20] as its impurity solver.
We summarize our two main results. First, the compressibility is clearly enhanced in the metallic phase close to the insulating phase, i.e., when the occupation number is slightly away from integer, near the critical temperature. In Ref. [5], such a compressibility enhancement (denoted as divergence therein) was first predicted for the one-and two-band Hubbard models (i.e., N = 2 and N = 4, respectively) near half filling and explained in terms of the Landau functional. In Ref. [7], the compressibility enhancement is observed for the N = 3 case. Here we generalize the scenario and provide a direct connection with spectral properties: a compressibility enhancement occurs near any integer occupation (except for completely empty and full occupation) for flavor-symmetric Hubbard models for general N , As temperature grows from 0 to T c , the compressibility gets enhanced, while at the same time the local spectral function gets suppressed by finite temperature.
Second, the quasiparticle peak in the metallic phase close to the Mott transition is necessarily always strongly asymmetric: the peak widths on positive and negative energy sides are different. This asymmetry arises from the strong imbalance of doublon and holon occupations in the slightly doped regime. We substantiate this argument by studying the doublon and holon correlation functions as done in our previous studies [21,22] on the subpeaks at the inner edges of the Hubbard bands.
The rest of this paper is organized as follows. In Sec. II, we provide details on the Hamiltonian and the numerical methods. In Sec. III, we study the compressibility, the local spectral function, and the local correlation functions of doublons and holons. Sec. IV offers our conclusion.

A. SU(N ) Hubbard model
The SU(N ) Hubbard model is the Hubbard model with N flavors of fermions which are fully symmetric, that is, Coulomb interaction strength U , hopping amplitude t, and chemical potential are independent of particle flavor ν. Its Hamiltonian is given by where c iν annihilates a particle with flavor ν = 1, . . . , N at lattice site i,n i ≡ ν c † iν c iν is particle number operator, i, j means nearest-neighbour sites, and µ is the chemical potential with offset such that µ = 0 yields particle-hole symmetry. The Hamiltonian has U(1) charge and SU(N ) flavor symmetry.
The SU(N ) Hubbard models had been originally considered as effective descriptions of multi-band strongly correlated materials, where the SU(N ) flavor symmetry is an approximation. Recently, these models were realized in ultracold atom experiments [10,12], where the SU(N ) flavor symmetry is exact and N is tunable over the range 2 ≤ N ≤ 6.
We take the chemical potential µ in Eq. (1) to be uniform throughout the lattice. Results for uniform systems are useful also for studying inhomogeneous systems within the context of the local density approximation (LDA). A detailed analysis [23] shows that the LDA is a good approximation in studying both the occupation number profile (real-space distribution) and the time-offlight (momentum-space distribution) for the Hubbard model in a harmonic trap.

B. DMFT
Dynamical mean-field theory (DMFT) [17,18] has been successfully used to study the paramagnetic Mott transition. In the single-site setting of DMFT, the Hubbard model is mapped onto the single-impurity Anderson model (SIAM). There the impurity, representing a lattice site with Coulomb interaction, is coupled to a bath of non-interacting fermions, and the energy dependence of the impurity-bath hybridization function encodes correlation effects (e.g., a Mott gap) within the rest of the lattice. The self-consistent solutions of the SIAM describe homogeneous phases of the original lattice model.
The mapping onto the effective impurity model relies on the approximation that the self-energy is local, i.e., momentum-independent, and charge or magnetic ordering is suppressed by assuming a fully frustrated lattice. This approximation of locality becomes exact in the limit of infinite coordination number of lattice z → ∞ [24]. To have finite bandwidth in this limit, the hopping amplitude is scaled as t ∝ 1/ √ z. Then the Green's function in the lattice is derived from the impurity Green's function, using the density of states ρ 0 (ω) of non-interacting lattice. In this work, we consider the semi-elliptic choice ρ 0 (ω) = 2 where D ≡ 2t √ z := 1 is the half-bandwidth which we set as the unit of energy Note that the Fermi energy of the non-interacting system is D + µ. We also set = k B = 1 throughout.
Due to this mapping, the overall feasibility as well as the accessible parameter range of DMFT depend on which method is used as impurity solver to solve the effective impurity model. Here we use the NRG as impurity solver, since it can provide the correlation functions on the real-frequency axis directly, thus avoiding the numerically ill-posed problem of having to analytically continue imaginary-frequency data to the real axis. Also, NRG is applicable to arbitrary temperature, including zero temperature at comparable computational cost. See the next subsection for details of the NRG method.

C. NRG
We solve the effective SIAM by using the full-densitymatrix NRG [25,26]. The bath is discretized on a logarithmic energy grid set by the coarse-graining parameter Λ = 4. The resulting discrete impurity model is mapped exactly then onto a Wilson chain with exponentially decaying hopping. By using energy scale separation, the iterative diagonalization yields a complete basis of approximate many-body eigenstates [27,28]. Here we keep up to N keep = 2500 low-energy multiplets at each of the early iterations corresponding to large energy scales. In later iterative diagonalization steps in the strongcoupling fixed point regime, for computational efficiency, we also apply a rescaled truncation energy threshold of E trunc = 9, which is expected to give converged results with keeping less multiplets than N keep [26]. Using a complete basis of energy eigenstates, the correlation functions at the impurity are obtained in the Lehmann representation as a collection of discrete spectral weights. To recover continuous spectral functions, we broaden the discrete spectral data with appropriate broadening kernels [25,29].
To simulate the multi-flavor SIAM with feasible computational cost, we exploit the U (1) charge ⊗ SU(N ) flavor symmetry in the system by making use of the QSpace tensor library for general non-Abelian symmetries [26,30]. This organizes the Hilbert space in terms of SU(N ) multiplets, and operates systematically at the level of reduced matrix elements, with the Clebsch-Gordan coefficients split off and dealt with separately. It allows us to efficiently perform DMFT+NRG calculations on multi-flavor models with SU(N ) symmetry up to N = 5, bearing in mind that typical multiplet sizes grow exponentially in N [30]. Furthermore, we use the adaptive broadening scheme [29] to improve the spectral resolution of correlation functions at higher energy. Specifically, we average the results over two discretization grids (n z = 2), followed by an adaptive log-Gaussian broadening whose width σ is controlled by the overall prefactor α = 2 and a lower bound σ ≥ (ln Λ)/8. At or below the energy scale of temperature T a linear broadening is further applied to smooth out artifacts at |ω| T . See Ref. [29] and the Supplementary Material of Ref. [21] for details.
Since the NRG calculation requires less computational cost for larger Λ, here we choose a rather large value Λ = 4 to explore O(10 4 ) data points of (N, U, T, µ) efficiently. As a tradeoff, we have limited spectral resolution at finite energies in the local correlation functions which is only partly regained by using the so-called z-averaging procedure, standard for NRG applications [20,31], with n z = 2. In particular, with Λ = 4 the discretization is too crude, e.g., to resolve the subpeaks that are known to occur at the inner Hubbard band edges [21,22]. However, such features are irrelevant for this work: The occupa- And the physical phases are rather determined by the low-energy part of A(ω), e.g., the quasiparticle peak or the Mott gap. Moreover, we expect that the doublon-holon subpeaks [21,22] are fairly suppressed in the vicinity of the Mott transition because of the strong asymmetry in the doublon-holon dynamics to be discussed below.
In this work, we study the occupation number n ≡ n i , the compressibilityκ [cf. Eq. (2) below] and local correlation functions A XY (ω) ≡ − 1 π Im X iν ||Y iν ω , which are the imaginary part of retarded time correlators of local operators X and Y acting on site i, transformed to the frequency domain. Based on the DMFT mapping onto the SIAM and the semi-elliptic density of states ρ 0 (ω), these local properties at a lattice site are equivalent to the same properties at the impurity when the self-consistent solution of the SIAM is achieved.

D. Compressibility
The compressibility is defined as where we only use the derivativeκ, without rescaling, for the remainder of the paper [12]. We obtain n for a linear grid of µ with grid size ∆µ = 0.05, and computẽ κ by numerically differentiating n(µ). Since the latter is sensitive to numerical noise, even if the curves n(µ) look smooth except at phase transition points (cf. Fig. 1), we determine the slope of n(µ) at µ = µ by fitting at most five consecutive points on the short interval µ − 2∆µ ≤ µ ≤ µ +2∆µ with a quadratic polynomial. When U is larger than the critical strength, a Mott transition occurs, signaled by discontinuities in n(µ) and/orκ(n), as discussed in much detail in the following sections. Thus, when µ is close to these values, we exclude the points n(µ) beyond the critical value to keep the fitting error minimal.

A. Zero temperature
We start with studying the filling-driven Mott transition at T = 0, which by definition is free from thermal fluctuations. These results are obtained by directly solving the SIAM at infinitesimally low temperature T = 0 + , i.e., not by extrapolating finite-T results [6]. Indeed, this  , 7) is not visible, since the coexistence region [µc1, µc2] in the curve n(µ) is narrower than the grid size ∆µ = 0.05. To indicate the second-order nature of MIT at T = 0, we have extrapolated the shaded areas to integer n, while the lines with downward arrows with finite slopes connect actual data points. The curves for n > N/2 can be deduced by particle-hole symmetryκ(n) =κ(N − n). accessibility of low temperatures is a major strength of using NRG as the DMFT impurity solver [33].
In Fig. 1 we present our data on the particle number per site n vs. the full range of chemical potential µ, and in Fig. 2 we depict the corresponding compressibilityκ vs. filling n/N . For weak interaction U = 1, the compressibilityκ is almost independent of n for n/N 0.1. In contrast,κ(n) has local minima at integer n for larger U ≥ 3. In Fig. 1, lines n(µ) with finite slope (i.e.,κ > 0) correspond to metallic phases. Conversely, horizontal plateaus represent the incompressible phase of a Mott insulator. These plateaus appear due to an interactiondriven Mott transition, i.e., by increasing U beyond a critical interaction strength U c (N, [n]) that depends on both the number of flavors N [34][35][36] as well as the integer filling [n] of the Mott plateau [4].
At each end of a plateau in the n(µ) curves (except for n = 0, N ), a Mott transition occurs, accompanied with a hysteresis loop [5,6], as shown (in the insets) in Fig. 1(c). For each hysteresis loop, we can define a pair of critical values µ c1 and µ c2 of the chemical potential: µ c1 is the value of the chemical potential at the outer edge of a plateau in n(µ) which thus describes an insulator-tometal transition (IMT). Similarly, µ c2 is the value where the metallic solution terminates within a Mott plateau, and thus describes a metal-to-insulator transition (MIT). Therefore, in between two critical values (µ c1 < µ < µ c2 and µ c2 < µ < µ c1 for the left and right ends of the Mott plateau, respectively), both insulating and metallic phases coexist, i.e., n(µ) is double-valued.
The compressibilityκ vs. n also shows discontinuities at integer n, associated with the Mott transition, as can be observed in Fig. 2(b) and (c). An IMT, depicted by an upward arrow, involve not only a jump inκ but also in the occupation n. Similarly, also across a MIT, depicted by a downward arrow, a jump in bothκ and n can occur. However, at T = 0, the jump in n should disappear, such that n evolves continuously across MIT [6]. Thus the downward arrows in Fig. 2 should in principle be strictly vertical; the reason why they are slanted, instead, is the nonzero grid size, ∆µ = 0.05, used for our calculations.
The continuity of n and the discontinuity ofκ across the metal-to-insulator transition reflect the second-order nature of the Mott transition at T = 0. Within the coexistence regime, there exists another critical value of chemical potential µ c at which the metallic and insulating solutions have the same values of free energy. For finite temperature 0 < T < T c , the transition at µ = µ c is first-order. In contrast, for T = 0, one has µ c = µ c2 = µ c1 and the transition at µ = µ c turns into a secondorder transition [5,6], in that n is continuous butκ is discontinuous (note that n andκ are proportional to the first and second derivatives of the free energy with respect to chemical potential [6]). Another exceptional situation for which the transition at µ = µ c is not first-order arises at the critical end point T = T c , where the coexistence region shrinks to a point, i.e., µ c = µ c1 = µ c2 , and, again, the transition becomes second-order.
For n = 0 and n = N the system is a band insulator, and thus no longer a Mott insulator. Correspondingly, we also observe no phase coexistence near the plateaus for n = 0 and N . Therefore while these trivial phases are still incompressible, their plateaus are excluded from our discussion of Mott plateaus. The value for the chemical potential below which the system becomes empty, is given by µ ≤ µ 0 ≡ −(N − 1) U 2 − D, in agreement with the overall trend seen in Fig. 1. This value can be motivated as follows: By substituting µ = −(N − 1) U 2 to µ, the first term in Eq. (1) favors zero and one occupation numbers equally. By adding another shift in the chemical potential, µ = −D, when associated with the second term in Eq. (1), this empties this non-interacting kinetic part of the Hamiltonian. Therefore the system becomes empty (n = 0) for µ ≤ µ 0 ≡ µ + µ , resulting in the above expression. Similarly, by particle-hole transformation, the system becomes completely filled (n = N ) for µ ≥ −µ 0 . appearing as discontinuous jumps inκ. Thick horizontal bars above the curves κ(n) (solid lines) near integer n (vertical gray lines) indicate the range of n for the metallic phase in the coexistence regime, color matched with the curvesκ for the same T . As T increases from 0 to Tc, the coexistence regions get narrower in n. For some curves, e.g., for N = 2, the coexistence region is narrower than the numerical grid size ∆µ = 0.05. As T is increased towards Tc, local maxima ofκ(n) appear near integer n, which we call compressibility enhancement. At T > Tc, the discontinuity in the curvesκ(n) disappears, which indicates that a crossover (rather than a phase transition) occurs between metallic and insulating behavior.

B. Finite T and compressibility enhancement
Next we analyze the effect of finite temperature T on the compressibility curvesκ(n), as shown in Fig. 3. We observe that the jumps inκ(n) near integer n survive for temperatures below a critical value T c , and disappear above it. That is, T c is the critical temperature of the paramagnetic Mott transition, and it depends on N , U , and µ. For the U = 7 considered in Fig. 3, we find that T c lies between 0.01 and 0.03.
Near and across the critical temperature, T ∼ T c , we observe a compressibility enhancement:κ(n) exhibits lo- cal maxima for n close, but not equal to, integer values. These local maxima ofκ(n) become more pronounced as T gets closer to T c , both from above and below. For example, for the curves of N = 5 and T ∈ [0.01, 0.03], a peak ofκ(n) associated with the compressibility enhancement is even the global maximum, not only a local maximum. In contrast, for T = 0 and T = 0.1, which are far below and above T c , respectively, the curvesκ(n) decrease monotonically as n approaches an integer both from above or below, and reach zero either by a jump for T < T c or continuously for T > T c .
The compressibility enhancement directly originates, by definition, from qualitative changes in the curves n(µ) for different T . In Fig. 4, we plot n(µ) for the same choice of parameter sets (N, U, T ) as in Fig. 3, but zooming in towards the coexistence region of the Mott transitions, choosing n slightly larger than N/2 , i.e., the downrounded value of N/2. As T increases from 0 to T c , the inner edge of the coexistence region, µ c2 , rapidly shifts towards the outer edge, µ c1 , while µ c1 likewise shifts outward, but slower than µ c2 does. Hence the width of the coexistence regime, |µ c1 − µ c2 |, decreases. In a sense, therefore, finite temperature destabilizes the metallicity, but stabilizes the insulting phase. Interestingly, this be- havior is similar to the T -dependence of U c2 and U c1 at half filling [14,36]. At the same time, the slope of the slanted part of n(µ) (corresponding to the metallic phase) in the close vicinity to integer n increases, leading to the compressibility enhancement in Fig. 3. We also see that contrary to the T = 0 case, not onlyκ but also n exhibits a jump at µ = µ c2 for 0 < T < T c .
Finally let us emphasize that while the data in Fig. 4 focusses on the particle-doped case (with n slightly above N/2 ), the overall behavior for the hole-doped case (n slightly below N/2 ) is completely analogous, bearing in mind that, via particle-hole transformation, one has n(µ) = N − n(−µ).

C. Local spectral functions
The coexistence region analyzed in Fig. 4 is intricately linked to a competition between the metallicity that permits non-integer occupation, and the Mottness that constrains n to be integer. Therefore in order to gain a better understanding, we now look into the local spectral functions in the metallic regime, with a focus on the mutual interplay between average local occupation n and the quasiparticle peak in the spectral data.
The decrease of the average local occupation n towards the Mott plateau [n] (the rounded value of n) as T is increased, as shown in Fig. 4, is necessarily connected to the thermal suppression of the quasiparticle peak in the local spectral function A(ω). Since its total weight, ∞ −∞ dω A(ω) = 1, is preserved by a sum rule, the ensuing transfer of spectral weight necessarily also influences the local occupation, n = N ∞ −∞ dω A(ω)/(e ω/T + 1). A detailed analysis of the spectral data in the metallic phase is presented in Fig. 5. There we show in each panel the local spectral functions for fixed N and µ but for several values of T . Because of the large U = 7, there are two well-developed Hubbard bands, centered at ω ∼ D − U and ω ∼ D, respectively. Since µ is chosen as a fixed value slightly larger than the critical values µ c2 for different T , the spectral data is strongly asymmetric around ω = 0 despite, e.g., n ≈ N/2 for even N . Specifically, the lower Hubbard band is well separated towards negative frequencies, whereas the lower edge of the upper Hubbard band is close to the Fermi level.
In addition to the Hubbard bands, the spectral functions in Fig. 5 for T ≤ 0.01 < T c feature a quasiparticle peak at the Fermi level ω = 0. As T increases, the quasiparticle peak gets suppressed and the occupation number n approaches [n] (see middle insets). The quasiparticle peaks for T < T c represent Fermi-liquid quasiparticles which, due to the narrowness of the peak, have heavy effective mass. This Fermi-liquid state hosts low-energy charge fluctuations, so non-integer n is generally possible [21].
On the other hand, for T ≥ 0.03 > T c , the significant suppression or the absence of the quasiparticle peak rules out coherent low-energy quasiparticles. So the state of the system is well described by the Hubbard bands only. If the lower Hubbard band (LHB) below the Fermi level is fully occupied and the upper one (UHB) above the Fermi level is empty, the lattice sites are filled by an integer number of particles, [n], without charge fluctuations. Accordingly, the integrated spectral weights of the individual Hubbard bands would be also integers, e.g., n LHB ≡ N LHB dω A(ω) = [n]. In Fig. 5, however where µ has been chosen slightly above µ c1 , the lower edge of the upper Hubbard band has dropped slightly below the Fermi level, thus making a small contribution to the occupancy. As a consequence, the values n(T = 0.03) are very close to, but slightly larger than, integers [n]. For even larger temperature T = 0.1, the thermal window of the Fermi-Dirac distribution (e ω/T + 1) −1 widens, so the occupation n deviates even more strongly from the integer [n].
The T -dependence of n and n LHB , presented in the middle insets of Fig. 5, show how the quasiparticle weight, i.e., the spectral weight of the quasiparticle peak, is transferred to the Hubbard bands as T increases. Since the total spectral weight is preserved, the difference n(0) − n LHB T =0 is equivalent to the negative-frequency part of the quasiparticle weight, i.e., integrated up to ω = 0. This weight is fully transferred to the Hubbard side bands as T is increased towards T c . This spectral weight transfer can be split into two net flows: the weight [n] − n LHB T =0 is transferred to the LHB, whereas the rest, i.e., n T =0 − [n], together with the quasiparticle weight for ω > 0, flows into the UHB. Surprisingly, despite the distance between the LHB and the Fermi level, a significant portion of the negative-frequency quasiparticle weight flows into the LHB. [e.g., see the change of height in the LHB in Fig. 5(d)].
The thermal suppression of the quasiparticle peak, which is accompanied by a transfer of spectral weight, thus pushes n closer towards [n] with increasing T at fixed µ. Correspondingly, µ c2 (T ) changes much more sensitively with T than µ c1 (T ). In the metallic phase near µ c2 (T ), the transfer of spectral weight from the quasiparticle peak into the LHB with increasing T necessarily leads to an increase in µ c2 , whereas no such weighttransfer occurs near µ c1 (T ) in the insulating phase, which lacks a quasiparticle peak. This is consistent with the fact, discussed in Sec. III B, that increasing T causes a stronger increase in µ c2 (T ) than in µ c1 (T ).
Once the chemical potential is outside the range of the Mott plateau (including the coexistence region), the system is always metallic, and the spectral functions evolve smoothly in terms of a crossover as temperature increases. Therefore in this case, quasiparticle peaks are present also for temperatures above T c , which denotes the critical temperature right at the IMT and MIT. This is the reason why quasiparticle peaks occur for T = 0.03 > T c in Fig. 5(a)-(b).
However, for these high-temperature peaks the evolution with decreasing µ at fixed T is qualitatively different from those of the low-temperature peaks in the regime T ≤ 0.01 < T c within the metallic phase. This is illustrated in Fig. 6, which shows how the local spectral function A(ω) evolves when µ is decreased (top to bottom) towards the MIT at µ c2 . As the MIT is approached while lowering µ at a given, fixed temperature, the quasiparticle peak behaves differently depending on whether that temperature lies above or below T c . In the former case, i.e., if the fixed temperature satisfies T ≥ 0.1 > T c , both the height and width of the quasiparticle peak smoothly decrease with decreasing µ, which is consistent with the crossover behaviour above T c . In contrast if the fixed temparature lies in the range T ≤ 0.01 < T c , the quasiparticle peak becomes narrower as µ decreases towards µ c2 , while its height remains almost unchanged. Once µ has passed below the MIT at µ c2 , the height A(ω = 0) drops abruptly, which is consistent with the transition nature below T c . In the limiting case of T = 0, the spectral function A(ω = 0) is pinned to the value 2/π all along the metallic phase, as dictated by the Luttinger theorem [37]. Our DMFT+NRG result in Fig. 5 fulfills this relation with accuracy better than 3% due to the intrinsic high accuracy of NRG at low energies, despite strongly broken particle-hole symmetry. For example, at T = 0 the curves in Fig. 5 have the zero-frequency values π 2 A(ω = 0) 0.9989, 0.9726, 0.9998, 0.9925 for N = 2, 3, 4, 5, respectively.
Again let us emphasize that, while the spectral data in Fig. 5 above is for the particle-doped case, the spectral functions for slight hole doping can be simply deduced by particle-hole transformation, which yields the equivalence A(ω)| µ = A(−ω)| −µ .
We briefly discuss the effect of a non-uniform potential. In Ref. [23], the paramagnetic Mott transition of the SU(2) Hubbard model has been studied by using real-space DMFT, which incorporates the non-uniformity of harmonic confinement potential. There the metallic phase was found to exist in a wider region than predicted by the LDA (which we use in this work; see Sec. II A), since the metallicity can "penetrate" into nearby insulating regions via the Kondo effect. On the other hand, the relation between the deviation of the local occupation number from integer and the thermal suppression of the quasiparticle peak was also found there, (e.g., see Fig. 5 of Ref. [23]), consistent with our result.

D. Doublon and holon correlators
The overall spectral data at finite doping as in Fig. 5, by construction, is always strongly particle-hole asymmetric. At weak doping and large U , the quasiparticle peak is necessarily close to one Hubbard band but clearly separated from the other. This asymmetry is also reflected in the shape of the quasiparticle peak itself.
It is possible to understand the origin of this asymmetry of the quasiparticle peaks by studying the correlation functions of doublons and holons [21]. We define a doublon (holon) as a local excitation that creates (annihilates) a particle at a lattice site filled by [n] particles. Accordingly the creation operators for doublons and holons are expressed as the projected operators (so-called Hubbard operators), where P i,[n] is the projector onto the subspace in which site i has [n] particles. For the N = 2 case, the fermion operator c iν can be decomposed into doublon and holon operators, c iν = d iν + h † iν . Thus the local spectral function can be decomposed as where the relation A h † d † (ω) = A dh (ω) holds generally. At the particle-hole symmetric point µ = 0, we have the further relation A dd † (ω) = A h † h (−ω). As before, the doublon and holon correlators for slightly hole doped cases can be deduced from the particle doping result in Fig. 7 via particle-hole symmetry: For N > 2, the decomposition of c iν involves other projected operators, such as c iν P i,[n]−1 , in addition to d iν and h † iν . In this case, the decomposition of A(ω) in Eq. (4a) is not exact but approximate, since additional terms arise beyond those shown in Eq. (4a). However, these additional terms are negligible for large U . Indeed, the deviation |A(ω) − X,Y =d,h † A XY † (ω)| becomes smaller, since the probability that a site contains less than [n]−1 or more than [n]+1 particles is suppressed due to the large cost in Coulomb energy. Therefore the projected particle operators which involve the projectors P i,m with m > [n]+1 or m < [n]−1 (i.e., other than doublon and holon operators) have negligible contribution to the correlation functions. Before analyzing specific results, let us discuss a few general properties of the correlators of doublons and holons. From Eq. (4a), we have three independent doublon and holon correlators, A dd † (ω), A dh (ω) = A h † d † (ω), and A h † h (ω), as shown in Fig. 7. They reflect the overall structure of the full spectral function, including the Hubbard bands and the quasiparticle peak in the metallic phase [21]. The doublon-doublon correlator A dd † and the holon-holon correlator A h † h have two major features. A dd † (A h † h ) carries most of the UHB (LHB), as well as the negative (positive) frequency side of the quasiparticle peak, respectively. The latter are centered around a small energy scale −ω s (ω s ), where ω s corresponds to the width of the quasiparticle peak centered around ω = 0, which itself can be related to the energy scale of flavorlike collective modes [22] via the (dynamical) flavor susceptibility. Note that, for N = 2, flavors equivalently represent spins for one band of electrons.
The features in the doublon-holon correlation functions necessarily correspond to dynamics at different energy scales associated with the Hubbard bands and the quasiparticle peaks, respectively. For simplicity, consider the case N = 2 at integer-filling n = 1 (the generalization for different N and [n] is straightforward). At T = 0, the positive and negative frequency sides of a correlator A XY directly correspond to the Fourier transforms of X(t)Y (0) and X † (t)Y † (0) * , respectively. Hence, for example, the contribution to the UHB by A dd † corresponds to the dynamics of a doublon excitation added at time 0 and then again removed at time t. Conversely, the low-energy feature centered at ω = −ω s < 0 means that a single spin remaining after removing a doublon at time 0 undergoes a time evolution governed by the spinlike collective mode with energy scale ω s , until a doublon is regenerated on top of the spin at time t. The features of A h † h can be explained analogously by swapping the roles of doublon and holon. On the other hand, the doublonholon correlators A dh (ω) = A h † d † (ω) mainly contribute to the quasiparticle peak, rather than to the Hubbard bands. It means that the doublon and holon excitations are combined at low energies to build quasiparticles. Now we demonstrate that the asymmetry of the quasiparticle peak in Fig. 5 originates from striking differences between the doublon and holon spectra, shown in Fig. 7. These differences stem from the strong asymmetry in energy cost for doublon and holon excitations, due to large U , despite the low level of particle doping (n−[n])/N < 0.03. In the slightly particle-doped regime, the UHB (LHB) originating from local doublon (holon) excitation lies close to (far from) the Fermi level. Due to this strong asymmetry, the metallic ground state in this particle-doped Mott insulator contains much more doublons than holons; correspondingly, A dd † (ω = 0) is higher than A h † h (ω = 0) by more than an order of magnitude. (This fact is evoked when constructing the t-J model [38], a widely used effective low-energy model for a Mott insulator, which in the case of hole-doping completely neglects double-occupancy.) In any case, the metallic ground state should contain delocalized doublon and/or holon fluctuations since otherwise the system would be insulating. Such fluctuations can optimize the total energy, by decreasing kinetic energy more than the related increase in Coulomb interaction energy. Finally, the differences in spectral strength between A dd † and A h † h in the quasiparticle peak regime, combined with the fact that the contributions of A dd † (A h † h ) to the quasiparticle peak are centered at ω = −ω s (+ω s ), necessarily cause the asymmetry of the quasiparticle peak in A(ω) [cf. Eq. (4a)].

IV. CONCLUSION
We have investigated the compressibility in the metallic and paramagnetic insulating phases along the fillingdriven Mott transition of the SU(N ) Hubbard model.