## 2 Numerical Methods and Preparations  
### 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:  
\[\hat{\mathcal{H}}=-t\sum_{\sigma=\uparrow,\downarrow}\sum_{\langle i,j\rangle} \left(\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{j,\sigma}+\text{h.c.}\right)+U\sum _{i}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}+\sum_{i,\sigma}\left[V(x_{i},y _{i})-\mu_{\sigma}\right]\hat{n}_{i,\sigma}\,, \tag{1}\]  
where \(\hat{c}^{\dagger}_{i,\sigma}\) and \(\hat{c}_{i,\sigma}\) are the creation and annihilation operators of the fermionic particle with a pseudospin component \(\sigma\) at site \(i\), respectively, and \(\hat{n}_{i,\sigma}\) is the corresponding number operator. The summation with \(\langle i,j\rangle\) runs over the nearest neighbor pairs of the sites \(i\) and \(j\). The physical constants \(\hbar\) 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)=\frac{1}{2}V_{0}(x^{2}+y^{2})\) with isotropic trap strength \(V_{0}\). The particle number \(N=N_{\uparrow}+N_{\downarrow}=\sum_{i,\sigma}\langle\hat{n}_{i,\sigma}\rangle\) and the total polarization \(P=(N_{\uparrow}-N_{\downarrow})/N\) can be controlled by tuning the chemical potential \(\mu_{\sigma}\). 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 \(\mathbf{r}\) in the trap is approximated from the solution of the translation invariant system by applying local chemical potential \(\mu^{\prime}=\mu-V(\mathbf{r})\), which is shifted by the trap potential at \(\mathbf{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 [74; 13; 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.