## 2 Numerical Methods and Preparations
### Real-Space Dynamical Mean-Field Theory
Let us begin with the diagonal matrix of the local self-energies \([\mathbf{\Sigma}(\sigma;\mathrm{i}\omega_{n})]_{ij}=\delta_{ij}\Sigma_{j}( \sigma;\mathrm{i}\omega_{n})\), where \(\omega_{n}\) is the fermionic Matsubara frequency, and \(\delta_{ij}\) is the Kronecker delta. The site-dependent \(\Sigma_{j}(\sigma;\mathrm{i}\omega_{n})\)is an initial guess at the first iteration. Then, given the self-energy matrix \(\mathbf{\Sigma}(\sigma;\mathrm{i}\omega_{n})\), the Dyson equation for the lattice Green's function \(\mathbf{G}(\sigma;\mathrm{i}\omega)\) can be written at each Matsubara frequency \(\omega_{n}\) as  
\[\mathbf{G}^{-1}(\sigma;\mathrm{i}\omega_{n})=\mathbf{G}_{0}^{-1}(\sigma; \mathrm{i}\omega_{n})-\mathbf{\Sigma}(\sigma;\mathrm{i}\omega_{n})\,, \tag{2}\]  
where \(\mathbf{G}_{0}\) is the lattice Green's function of the non-interacting system, calculated as  
\[[\mathbf{G}_{0}^{-1}(\sigma;\mathrm{i}\omega_{n})]_{ij}=[\mu_{\sigma}-V(x_{i},y_{i})+\mathrm{i}\omega_{n}]\delta_{ij}-[\mathbf{H}_{0}]_{ij}\,, \tag{3}\]  
for the trap potential \(V(x_{i},x_{j})\) and the hopping part \([\mathbf{H}_{0}]_{ij}\) between sites \(i\) and \(j\) in the original Hubbard Hamiltonian. Numerically solving Equation (2) to obtain \(\mathbf{G}(\sigma;\mathrm{i}\omega_{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\times 2700\) matrix of \(\mathbf{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 \(\mathbf{G}\) calculated from the lattice Dyson equation, one can generate the site-dependent Weiss mean-field from the local Dyson equation as:  
\[1/\mathcal{G}_{i}^{(0)}(\sigma;\mathrm{i}\omega_{n})=1/\left[\mathbf{G}(\sigma ;\mathrm{i}\omega_{n})\right]_{ii}+\Sigma_{i}(\sigma;\mathrm{i}\omega_{n})\,. \tag{4}\]  
The Anderson impurity model is constructed at each site by using \(\mathcal{G}_{i}^{(0)}\) 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:  
\[\mathcal{H}_{i,n_{s}}=\mathrm{U}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}+ \sum_{\sigma=\uparrow,\downarrow}[V(x_{i},y_{i})-\mu_{\sigma}]\hat{n}_{i, \sigma}+\sum_{\sigma=\uparrow,\downarrow}\sum_{\hat{a}=1}^{n_{s}}\left[E_{ \alpha,\sigma}\hat{a}_{\alpha,\sigma}^{\dagger}\hat{a}_{\alpha,\sigma}+V_{ \alpha,\sigma}\hat{a}_{\alpha,\sigma}^{\dagger}\hat{c}_{i,\sigma}+V_{\alpha, \sigma}^{*}\hat{c}_{i,\sigma}^{\dagger}\hat{a}_{\alpha,\sigma}\right]\,, \tag{5}\]  
where \(\hat{a}^{\dagger}\) and \(\hat{a}\) are the creation and annihilation operators of the bath orbitals, and the parameters \(E_{\alpha,\sigma}\) and \(V_{\alpha,\sigma}\) represent the energy of the bath orbitals and the coupling to the impurity site.