Qubit parametrization of the variational discrete action theory for the multiorbital Hubbard model

TL;DR

This paper develops a qubit parametrization of the variational discrete action theory (VDAT) at =3 for the multiorbital Hubbard model, enabling efficient ground state calculations and analytical insights into Mott and Hund physics.

Contribution

It introduces a qubit-based variational approach for VDAT at =3, improving computational efficiency and analytical understanding of complex multiorbital systems.

Findings

Successfully applied to models with 2-7 orbitals.

Derived analytical expressions for critical U in large N_{orb} limit.

Demonstrated equivalence to slave spin mean-field theory for certain variants.

Abstract

The variational discrete action theory (VDAT) at \mathcal{N}=3 is a potent tool for accurately capturing Mott and Hund physics at zero temperature in d=\infty at a cost comparable to the Gutzwiller approximation, which is recovered by VDAT at \mathcal{N}=2. Here we develop a qubit parametrization of the gauge constrained algorithm of VDAT at \mathcal{N}=3 for the multiorbital Hubbard model with general density-density interactions. The qubit parametrization yields an explicit variational trial energy, and the variational parameters consist of the momentum density distribution, the shape of a reference fermi surface, and the pure state of a qubit system with dimension of the local Hilbert space. To illustrate the power of the qubit parametrization, we solve for the ground state properties of the multiorbital Hubbard model with Hund coupling for local orbital number N_{orb}=2-7. A Taylor series expansion of the partially optimized trial energy is used to explain how the Hund's coupling changes the order of the Mott transition. For the case of the SU(2N_{orb}) Hubbard model, an explicit approach for computing the critical U_{c} for the Mott transition is provided, yielding an analytical expression for U_{c} in the large N_{orb} limit. Additionally, we provide an analytical solution for the ground state properties of the single band Hubbard model with a special density of states. Finally, we demonstrate that the qubit parametrization can also be applied to \mathcal{N}=2, for both G-type and B-type variants, where the G-type yields an identical expression to the slave spin mean-field theory. The qubit parametrization not only improves the efficiency and transparency of VDAT at \mathcal{N}=3, but also provides the key advances for the construction of a one-body reduced density matrix functional capable of capturing Mott and Hund physics.

Tables

Table 1. Table 1: The values of Δ \Delta , U U , and d d at various critical points for the two peak Hubbard model solved using ρ ^ G ​ 2 \hat{\rho}_{G2} , ρ ^ B ​ 2 \hat{\rho}_{B2} , and ρ ^ G ​ 3 \hat{\rho}_{G3} . For ρ ^ B ​ 2 \hat{\rho}_{B2} , the critical point for local stability ( ρ ^ B ​ 2 \hat{\rho}_{B2} (L)) and the transition to the Hartree-Fock solution ( ρ ^ B ​ 2 \hat{\rho}_{B2} (HF)) are provided. For ρ ^ G ​ 2 \hat{\rho}_{G2} and ρ ^ G ​ 3 \hat{\rho}_{G3} , the critical values for the Mott transition are provided.

	$\Delta$	$U/\left|K_0\right|$	$d$
${\hat{\rho }}_{G3}$	$\begin{array}{c}\frac{1}{6}\left(3-\sqrt{3}\right)\\ \approx 0.211325\end{array}$	$\begin{array}{c}\frac{3}{8}\left(9+5\sqrt{3}\right)\\ \approx 6.6226\end{array}$	$\begin{array}{c}\frac{1}{36}\left(32\sqrt{3}-55\right)\\ \approx 0.0118229\end{array}$
${\hat{\rho }}_{G2}$	$\frac{1}{4}$	$8$	$0$
${\hat{\rho }}_{B2}$ (L)	$\begin{array}{c}\frac{1}{12}\left(3-\sqrt{3}\right)\\ \approx 0.105662\end{array}$	$\frac{3\sqrt{3}}{2}\approx 2.59808$	$\frac{5}{36}\approx 0.138889$
${\hat{\rho }}_{B2}$(HF)	$\frac{1}{6}$	$\frac{27}{8}=3.375$	$\frac{17}{324}\approx 0.0524691$

Table 2. Table 2: The values of Δ \Delta , U U , and d d at various critical points for the single orbital Hubbard model on the Bethe lattice in d = ∞ d=\infty solved using ρ ^ G ​ 2 \hat{\rho}_{G2} , ρ ^ B ​ 2 \hat{\rho}_{B2} , and ρ ^ G ​ 3 \hat{\rho}_{G3} . For ρ ^ B ​ 2 \hat{\rho}_{B2} , the critical point for local stability ( ρ ^ B ​ 2 \hat{\rho}_{B2} (L)) and the transition to the Hartree-Fock solution ( ρ ^ B ​ 2 \hat{\rho}_{B2} (HF)) are provided. For ρ ^ G ​ 2 \hat{\rho}_{G2} and ρ ^ G ​ 3 \hat{\rho}_{G3} , the critical values for the Mott transition are provided.

	$\Delta$	$U/\left|K_0\right|$	$U/t$	$d$	$A$
${\hat{\rho }}_{G3}$	$\begin{array}{c}\frac{1}{6}\left(3-\sqrt{3}\right)\\ \approx 0.211325\end{array}$	$6.61836$	$5.61784$	$0.0162523$	$0.491668$
${\hat{\rho }}_{G2}$	$\frac{1}{4}$	$8$	$6.79061$	$0$	$\frac{1}{2}$
${\hat{\rho }}_{B2}$ (L)	$0.107273$	$2.60206$	$2.2087$	$0.166367$	$0.380259$
${\hat{\rho }}_{B2}$ (HF)	$0.153687$	$3.04071$	$2.58103$	$0.089871$	$0.447304$

Table 3. Table 3: Critical values for the SU(2 N o ​ r ​ b N_{orb} ) Hubbard model ( J / U = 0 J/U=0 ) on the Bethe lattice in d = ∞ d=\infty solved using ρ ^ G ​ 3 \hat{\rho}_{G3} .

$N_{orb}$	${\Delta }_c$	$b_c/t$	$A_c$	$U_c/t$
$1$	$\frac{1}{6}\left(3-\sqrt{3}\right)\approx 0.211325$	$5.34164$	$0.491668$	$5.61784$
$2$	$\frac{1}{20}\left(10-\sqrt{30}\right)\approx 0.226139$	$8.80352$	$0.496835$	$8.97283$
$3$	$\frac{1}{2}-\frac{1}{\sqrt{14}}\approx 0.232739$	$12.2288$	$0.498345$	$12.3511$
$4$	$\frac{1}{12}\left(6-\sqrt{10}\right)\approx 0.236477$	$15.6412$	$0.498984$	$15.7369$
$5$	$\frac{1}{22}\left(11-\sqrt{33}\right)\approx 0.238884$	$19.0475$	$0.499314$	$19.1261$
$6$	$\frac{1}{2}-\frac{\sqrt{\frac{7}{26}}}{2}\approx 0.240563$	$22.4505$	$0.499505$	$22.5172$
$7$	$\frac{1}{2}-\frac{1}{\sqrt{15}}\approx 0.241801$	$25.8515$	$0.499627$	$25.9094$

Table 4. Table 4: ℱ ​ ( n i , n j ) \mathcal{F}\left(n_{i},n_{j}\right) and E i ​ n ​ t E_{int} for the HF, MBB, Power, CA, and CGA functionals for the single orbital Hubbard model at half-filling in the paramagnetic phase.

Functional	$ℱ(n_i,n_j)$
HF	$n_i{n}_j$
MBB	$\sqrt{n_i{n}_j}$
power	${\left(n_i{n}_j\right)}^{\alpha }$
CA	$n_i{n}_j+\sqrt{n_i\left(1-n_i\right)n_j\left(1-n_j\right)}$
CGA	$\frac{1}{2}\left(n_i{n}_j+\sqrt{n_i\left(2-n_i\right)n_j\left(2-n_j\right)}\right)$
	$E_{int}$
HF	$\frac{U}{4}$
MBB	$\frac{U}{2}\left(1-2{\left(\int 𝑑k\sqrt{n_k}\right)}^2\right)$
power	$\frac{U}{2}\left(1-2{\left(\int 𝑑k{n}_k^{\alpha }\right)}^2\right)$
CA	$U\left(\frac{1}{4}-{\left(\int 𝑑k\sqrt{n_k\left(1-n_k\right)}\right)}^2\right)$
CGA	$U\left(\frac{3}{8}-\frac{1}{2}{\left(\int 𝑑k\sqrt{n_k\left(2-n_k\right)}\right)}^2\right)$

Equations

$$\hat{H}={\hat{H}}_K+{\hat{H}}_{loc}=\sum _{k\ell }{ϵ}_{k\ell }{\hat{n}}_{k\ell }+\sum _{i,\ell <{\ell }^{\prime }}{U}_{\ell {\ell }^{\prime }}{\hat{n}}_{i\ell }{\hat{n}}_{i{\ell }^{\prime }},$$

$$\hat{\rho }={\hat{K}}_1{\hat{P}}_1\dots {\hat{K}}_{\mathcal{N}}{\hat{P}}_{\mathcal{N}},$$

$${\hat{X}}_{i\Gamma }=\prod _{\ell =1}^{2N_{orb}}{\hat{X}}_{i\Gamma ;\ell },$$

$${\hat{X}}_{i\Gamma ;\ell }=\left(1-{\hat{n}}_{i\ell }\right){\delta }_{0,\Gamma \left(\ell \right)}+{\hat{n}}_{i\ell }{\delta }_{1,\Gamma \left(\ell \right)},$$

$$\displaystyle\hat{K}_{\tau}=\prod_{k\ell}\big{(}(1-\lambda_{k\ell\tau})(1-\hat{n}_{k\ell})+\lambda_{k\ell\tau}\hat{n}_{k\ell}\big{)},$$

$${\hat{P}}_{\tau }=\prod _i\left(\sum _{\Gamma }{u}_{\Gamma \tau }{\hat{X}}_{i\Gamma }\right),$$

$${\hat{\rho }}_{G2}={\hat{P}}_1{\hat{K}}_2{\hat{P}}_1^{†}.$$

$${\hat{\rho }}_{B2}={\hat{K}}_1{\hat{P}}_1{\hat{K}}_1^{†}.$$

$${\hat{\rho }}_{G3}={\hat{K}}_1{\hat{P}}_1{\hat{K}}_2{\hat{P}}_1^{†}{\hat{K}}_1^{†}.$$

$$e^{-\beta H}\approx {\left(e^{-\frac{\beta }{\mathcal{N}}\left({\sum }_{k\ell }{ϵ}_{k\ell }{\hat{n}}_{k\ell }\right)}{e}^{-\frac{\beta }{\mathcal{N}}\left({\sum }_{i,\ell <{\ell }^{\prime }}{U}_{\ell {\ell }^{\prime }}{\hat{n}}_{i\ell }{\hat{n}}_{i{\ell }^{\prime }}\right)}\right)}^{\mathcal{N}},$$

$${\mathbf{g}}^{-1}-1=\left({\mathbf{g}}_0^{-1}-1\right)\mathbf{S}$$

$${\mathbf{S}}_0={\mathbf{S}}_Q{\mathbf{S}}_K,$$

$${\mathbf{S}}_F={\mathbf{S}}_Q{\mathbf{S}}_K\mathbf{S},$$

$${\mathbf{g}}_0={\left(1+{\mathbf{S}}_0\right)}^{-1},$$

$$\mathbf{g}={\left(1+{\mathbf{S}}_F\right)}^{-1}.$$

$$e^{-{\mathbf{\mu }}_{\tau }\cdot \hat{\mathbf{n}}}{\hat{a}}_{\ell }^{†}{e}^{{\mathbf{\mu }}_{\tau }\cdot \hat{\mathbf{n}}}=\sum _{{\ell }^{\prime }}{\left[e^{-{\mathbf{\mu }}_{\tau }^T}\right]}_{\ell {\ell }^{\prime }}{\hat{a}}_{{\ell }^{\prime }}^{†},$$

$$e^{-{\mathbf{\mu }}_{\tau }\cdot \hat{\mathbf{n}}}{\hat{a}}_{\ell }{e}^{{\mathbf{\mu }}_{\tau }\cdot \hat{\mathbf{n}}}=\sum _{{\ell }^{\prime }}{\hat{a}}_{{\ell }^{\prime }}{\left[e^{{\mathbf{\mu }}_{\tau }^T}\right]}_{{\ell }^{\prime }\ell }.$$

$${\mathbf{g}}^{\prime }={\mathbf{N}}_b\mathbf{g}{\mathbf{N}}_b^{-1},$$

$${\mathbf{S}}_F^{\prime }={\left({\mathbf{g}}^{\prime }\right)}^{-1}-1={\mathbf{N}}_b{\mathbf{S}}_F{\mathbf{N}}_b^{-1},$$

$${\mathbf{S}}_K^{\prime }={\tilde{\mathbf{N}}}_b{\mathbf{S}}_K{\mathbf{N}}_a={\mathbf{S}}_Q^{-1}{\mathbf{N}}_b{\mathbf{S}}_Q{\mathbf{S}}_K{\mathbf{N}}_a.$$

$${\mathbf{S}}_0\mathbf{S}={\mathbf{S}}_Q{\mathbf{S}}_K\mathbf{S}={\mathbf{S}}_F,$$

$${{\mathbf{S}}^{\prime }}_0{\mathbf{S}}^{\prime }={\mathbf{S}}_Q{\mathbf{S}}_K^{\prime }{\mathbf{S}}^{\prime }={\mathbf{S}}_F^{\prime }.$$

$${{\mathbf{S}}^{\prime }}_0={\mathbf{S}}_Q{\mathbf{S}}_K^{\prime }={\mathbf{N}}_b{\mathbf{S}}_Q{\mathbf{S}}_K{\mathbf{N}}_a={\mathbf{N}}_b{\mathbf{S}}_0{\mathbf{N}}_a,$$

$${\mathbf{S}}^{\prime }$$

$$={\mathbf{N}}_a^{-1}\mathbf{S}{\mathbf{N}}_b^{-1},$$

$${\mathbf{g}}_0^{\prime }={\left(1+{\mathbf{S}}_0^{\prime }\right)}^{-1}={\left(1+{\mathbf{N}}_b\left({\mathbf{g}}_0^{-1}-1\right){\mathbf{N}}_a\right)}^{-1}.$$

$$\text{@text@baccentspdsymb^}={\text{@text@baccentspdsymb^}}_0\text{@text@baccentP^},$$

$$⟨{\text{@text@baccentH^}}^{(\mathcal{N})}{⟩}_{\text{@text@baccentspdsymb^}}=⟨{\text{@text@baccentH^}}_K^{\left(\mathcal{N}\right)}{⟩}_{{\text{@text@baccentρ^}}_K}+⟨{\text{@text@baccentH^}}_{loc}^{\left(\mathcal{N}\right)}{⟩}_{{\text{@text@baccentρ^}}_{loc}},$$

$${\text{@text@baccentρ^}}_K$$

$${\text{@text@baccentρ^}}_{loc}$$

$$\left({\mathbf{g}}_i^{-1}-1\right)=\left({\mathcal{G}}_i^{-1}-1\right){\mathbf{S}}_i,$$

$${\text{@text@baccentρ^}}_{loc;i}=\exp \left(-\ln {\left({\mathcal{G}}_i^{-1}-1\right)}^T\cdot {\text{@text@baccentn^}}_i\right){\text{@text@baccentP^}}_i.$$

$${\mathbf{g}}_i={\mathbf{g}}_i^{\prime },$$

$$⟨\text{@text@baccentO^}{⟩}_{{\text{@text@baccentρ^}}_{loc;i}}=\frac{⟨{\text{@text@baccentP^}}_i\text{@text@baccentO^}{⟩}_{{\text{@text@baccentρ^}}_{loc;i,0}}}{⟨{\text{@text@baccentP^}}_i{⟩}_{{\text{@text@baccentρ^}}_{loc;i,0}}}=\frac{{\left(\text{@text@baccentO^}\right)}_u\cdot \left({\otimes }_{\tau }{u}_{\tau }\right)}{{\left(\text{@text@baccent1^}\right)}_u\cdot \left({\otimes }_{\tau }{u}_{\tau }\right)},$$

$${\text{@text@baccentρ^}}_{loc;i,0}=\exp (-\ln {\left({\mathcal{G}}_i^{-1}-1\right)}^T\cdot {\text{@text@baccentn^}}_i),$$

Full text

I Introduction

The multiorbital Hubbard model is a minimal model of interacting electrons which accounts for key elements of typical strongly correlated electron systems Imada19981039 ; Fernandes2017014503 , and the solution in infinite dimensions provides a natural starting point for understanding the solution in finite dimensions Georges199613 ; Kotliar2006865 . The de facto standard for accurately solving the multiorbital Hubbard model in $d=\infty$ is the dynamical mean-field theory (DMFT) Georges199613 ; Kotliar200453 ; Kotliar2006865 ; Vollhardt20121 , which requires self-consistently solving an effective Anderson impurity model. Therefore, the proficiency of DMFT is determined by the state-of-the-art for quantum impurity solvers Gull2011349 , which are still extremely limited for the multiorbital case at zero temperature (see introduction of Ref. Cheng2022205129 ). Alternatively, the standard for efficiently solving the multi-orbital Hubbard model is the Gutzwiller approximation (GA), which semi-quantitatively describes the Fermi liquid phase and the metal-insulator transition, but only produces a crude description of the insulating phase Lu19945687 ; Okabe19972129 ; Bunemann19974011 ; Bunemann199829 ; Bunemann2007193104 ; Lanata2008155127 . A long standing goal has been to improve upon GA while not incurring the computational cost of DMFT. Fortunately, the variational discrete action theory (VDAT) provides an alternate route to the exact ground state properties in $d=\infty$ Cheng2021195138 ; Cheng2021206402 ; Cheng2022205129 ; Cheng2023035127 . The variational ansatz of VDAT is the sequential product density matrix (SPD), which is characterized by an integer $\mathcal{N}$ and as either G-type or B-type, where G-type $\mathcal{N}=1$ recovers the Hartree-Fock wave function, G-type $\mathcal{N}=2$ recovers the Gutzwiller wave function, and $\mathcal{N}\rightarrow\infty$ recovers the exact solution. A remarkable fact is that the G-type $\mathcal{N}=3$ SPD accurately captures Mott and Hund physics in $d=\infty$ Cheng2022205129 ; Cheng2023035127 , while maintaining a computational cost comparable to the Gutzwiller wave function.

An SPD with arbitrary $\mathcal{N}$ can be exactly evaluated in $d=\infty$ using the self-consistent canonical discrete action theory (SCDA) Cheng2021195138 ; Cheng2021206402 . There are currently three algorithms for executing VDAT within the SCDA, with the first two being completely general and the third one having some restrictions. The first approach is the most straightforward, requiring iterations over two steps Cheng2021206402 . In the first step, given the variational parameters, the SCDA equations can be self-consistently solved, yielding the energy for the corresponding SPD. In the second step, the variational parameters must be updated to minimize the energy. This straightforward approach was successfully executed on the single band Hubbard model in $d=\infty$. The second approach is the decoupled minimization algorithm Cheng2022205129 , where one simultaneously minimizes the variational parameters and updates the SCDA equations towards self-consistency. This decoupled minimization algorithm was executed for the G-type $\mathcal{N}=2-4$ SPD’s in the two orbital Hubbard model with the full rotationally invariant local interactions, demonstrating that $\mathcal{N}=3$ accurately describes Mott and Hund physics over all parameter space. The third approach is the gauge constrained algorithm Cheng2023035127 , which was originally applied to the case of a G-type $\mathcal{N}=2$ and $\mathcal{N}=3$ SPD with kinetic projectors that are diagonal in momentum space and interacting projectors that do not introduce off-diagonal terms at the single-particle level. The key feature of the gauge constrained algorithm is that the SCDA self-consistency can be automatically satisfied, greatly simplifying the task of minimizing over the variational parameters. Furthermore, the momentum density distribution is used to reparametrize the variational parameters in the kinetic projectors, and the trial energy can be straightforwardly optimized over the momentum density distribution. It is useful to note that the GA can only evaluate a G-type $\mathcal{N}=2$ SPD under a Gutzwiller gauge Cheng2023035127 . The SCDA has no such restrictions, which is formally appealing, but imposing certain gauge constraints allows the SCDA self-consistency to be automatically satisfied, while not reducing the variational capacity.

One shortcoming of the original gauge constrained algorithm for the G-type $\mathcal{N}=3$ SPD is that some of the variational parameters are unintuitive and there are two linear constraints on the momentum density distribution per spin-orbital Cheng2023035127 . In this paper, we introduce a new parametrization of the gauge constrained algorithm, referred to as the qubit parametrization, which is mathematically equivalent to the original gauge constrained algorithm and offers several important advantages. First, there is only one linear constraint per spin-orbital on the momentum density distribution, and therefore the number of variational parameters is reduced by one per spin-orbital. Second, the variational parameters are physically intuitive and facilitate a deeper understanding of how the SPD captures Mott and Hund physics. Therefore, the qubit parametrization achieves the long sought goal of resolving the shortcomings of the Gutzwiller approximation while maintaining the computational simplicity and physical appeal. The convenience of the qubit parametrization has allowed for the construction of a one-body reduced density matrix functional for the multi-orbital Hubbard model which exactly encapsulates the VDAT result at $\mathcal{N}=3$, and this is presented in a companion manuscript companion .

It is useful to compare to existing approaches in the literature which seek to achieve our same goal of efficiently and accurately solving the ground state properties of multiorbital Hubbard models. One approach is the ghost Gutzwiller approximation (gGA) Lanata2017195126 , which introduces the Gutzwiller wave function in an extended Hilbert space where the GA is then applied, though this comes with a substantial increase in computational cost over the usual GA given that ghost orbitals are introduced. In $d=\infty$, the gGA exactly evaluates the Gutzwiller wave function in the extended Hilbert space, and therefore provides an upper bound on the exact energy with increasing accuracy as the extended Hilbert space grows. The gGA has not been proven to converge to the exact solution of the multiorbital Hubbard model with increasing number of ghost orbitals, though numerical evidence suggests that this might be the case Lee2023L121104 ; Lee2023245147 . Empirically, it has been demonstrated that at least two ghost orbitals per physical orbital are needed for a global improvement over the GA Lanata2017195126 ; Lee2023245147 ; Lee2023L121104 , which substantially increases the computational cost of the gGA relative to GA. Therefore, the important practical comparison between VDAT and gGA contrasts VDAT at $\mathcal{N}=3$ to gGA with two ghost orbitals per physical orbital. VDAT at $\mathcal{N}=3$ has been straightforwardly applied up to the eight orbital Hubbard model Cheng2023035127 , while the gGA has been applied up to the three orbital model MejutoZaera2023235150 ; Lee2023245147 ; Giuli2025020401 using exact diagonalization and the five orbital model Lee2024115126 using density matrix renormalization group (DMRG). The success of the gGA in the multiorbital Hubbard model critically depends on the ability of DMRG or related techniques to approximately obtain a ground state in an expanded local Hilbert space. In terms of accuracy, both approaches appear to be excellent over all parameter space, though it seems possible that VDAT is more adept at capturing subtle differences between competing phases (e.g. compare Fig. 4 in Ref. Cheng2022205129 with Fig. 2 in Ref. MejutoZaera2023235150 ).

Another attempt to resolve the limitations of the GA is the slave boson approach Kotliar19861362 , which formally recasts the Hubbard model as a Hamiltonian of constrained fermions and auxiliary bosons. The simplest saddle point approximation is typically referred to as slave boson mean-field theory (SBMF), which exactly recovers the GA in the multi-orbital Hubbard model Bunemann2007193104 and has extensive applications Lechermann2007155102 ; Isidori2009115120 ; Bunemann2011203 ; medici2017167003 ; Piefke2018125154 ; Isidori2019186401 ; Chatzieleftheriou2020205127 . Given that the slave boson approach is an exact formalism, there have been various results including fluctuations beyond the mean-field solution Lavagna1990142 ; Jolicoeur19912403 ; Li1991369 ; Raimondi199311453 ; Arrigoni19933178 ; Li199417837 ; Zimmermann199710097 , though none of these studies address corrections to the ground state properties in the thermodynamic limit. A related idea is the time dependent Gutzwiller approximation Seibold20012605 ; Lorenzana2003066404 ; Schiro2010076401 ; Oelsen2011076402 ; oelsen2011113031 ; Schiro2011165105 ; Bunemann2015550 ; Fabrizio2017075156 , which has mainly focussed on computing response functions. One study did demonstrate that the ground state properties of the single orbital Hubbard model can be improved by updating the double occupancy based on the density-density response function Seibold20012605 , but this would be impractical in the multi-orbital Hubbard model, and we are not aware of any such applications. Therefore, neither slave bosons nor time dependent GA has provided a path beyond the GA for ground state properties, though the idea of ghost orbitals has recently been adopted into the slave boson methodology Lanata2022045111 , where the saddle point solution recovers the gGA. Another approach related to SBMF is the slave spin approach De'medici2005205124 , which recasts the Hubbard model with density-density interactions as a Hamiltonian of constrained fermions and auxiliary spins, and the saddle point approximation also recovers the GA. As in the case of slave bosons, the slave spin approach has produced many interesting results on the multi-orbital Hubbard model Hassan2010035106 ; medici2017167003 ; Georgescu2015235117 ; Maurya2021425603 ; Maurya2022055602 ; Crispino2023155149 ; Gorni2023125166 , but has not offered a path beyond the mean-field solution. In the present work, we demonstrate that the qubit parametrization for a G-type $\mathcal{N}=2$ SPD is identical to the slave spin mean field theory (see Appendix B).

Finally, it is useful to discuss the off-shell effective energy theory (OET) Cheng2020081105 , which can be reinterpreted through the G-type and B-type SPD for $\mathcal{N}=2$. OET introduces an approximation to evaluate the energy under the SPD, known as the central point expansion (CPE), and this is applied to both the G-type and B-type SPD for $\mathcal{N}=2$. Though the CPE is intrinsically different than the SCDA, it also exactly evaluates the $\mathcal{N}=2$ SPD in $d=\infty$, which is proven in the present work (see Appendix C). OET introduces corrections to the CPE for the G-type and B-type SPD, guaranteeing the limiting behaviors for the weak and strong coupling limits, respectively. Subsequently, one evaluates the total energy in both cases and selects the solution with lower energy. OET has yielded accurate results for the single band Hubbard model in $d=1$, $d=2$, and $d=\infty$. The correction to the CPE becomes far more challenging in the multiorbital Hubbard model, which has not been pursued further. However, a similar idea of empirically correcting the Gutzwiller approximation has been developed in the context of the correlation matrix renormalization approximation Yao2014045131 ; Liu2021081113 ; Liu2021095902 ; Liu2022205124 .

The paper is organized in the following manner. Section II provides an overview of the SPD ansatz, the gauge symmetry of the SPD, the SCDA, and the tensor product representation for evaluating expectation values. Additionally, several methodological generalizations are provided. Section III provides a high level overview of the qubit energy form for the G-type and B-type SPD at $\mathcal{N}=2$ and the G-type SPD at $\mathcal{N}=3$. Detailed derivations of the results from Section III are presented in Sections IV, V, and VI. Finally, the qubit energy form is used to examine the multi-orbital Hubbard model at half filling with $J/U\geq 0$ in Section VII.

II Review of the variational discrete

action theory

Here we review several key ingredients of VDAT, including the sequential product density matrix (SPD), the self-consistent canonical discrete action theory (SCDA), and the tensor product representation for evaluating local expectation values within the SCDA. In addition to reviewing these concepts, we generalize the gauge symmetry of the SPD beyond the case of a diagonal transformation, and generalize the tensor product representation to arbitrary $\mathcal{N}$. These developments will facilitate the subsequent analysis in this work.

II.1 The sequential product density matrix

The SPD is the variational ansatz used within VDAT Cheng2021206402 ; Cheng2021195138 , and the SPD is a straightforward generalization of a variational wave function ansatz first articulated in Ref. Dzierzawa19951993 . In this work, we focus on the multiorbital Hubbard model, given as

[TABLE]

where $k$ and $i$ denote momentum and real-space site indices, respectively, $\ell=1,\dots,2N_{orb}$ is the spin-orbital index within a local site, and $U_{\ell\ell^{\prime}}$ denotes the Coulomb interaction. In the context of the multiorbital Hubbard model, the SPD with an integer time step $\mathcal{N\geq}1$ is given as

[TABLE]

where $\hat{K}_{\tau}=\exp\left(\sum_{k\ell}\gamma_{k\ell\tau}\hat{n}_{k\ell}\right)$ is the kinetic projector, $\hat{P_{\tau}}=\exp(\sum_{i\Gamma}\upsilon_{\Gamma\tau}\hat{X}_{i\Gamma})$ is the interacting projector, $\tau=1,\dots,\mathcal{N}$ is the integer time index, the diagonal Hubbard operator $\hat{X}_{i\Gamma}$ is defined as

[TABLE]

and $\Gamma\left(\ell\right)$ is $\ell$-th bit in the binary representation of $\Gamma-1$, given as $\left(\Gamma\left(1\right)\dots\Gamma\left(2N_{orb}\right)\right)_{2}=\Gamma-1$, and $\Gamma=1,\dots,2^{2N_{orb}}$. It should be noted that in general, $\hat{P}_{\tau}$ is a direct product of local projectors from all sites while $\hat{K}_{\tau}$ is a general non-interacting projector. There are two schemes to ensure that the SPD is Hermitian and positive semi-definite, denoted as G-type or B-type Cheng2021195138 . In the following, we will examine both types for $\mathcal{N}\leq 3$. For simplicity, we restrict $\gamma_{k\ell\tau}$ and $\upsilon_{\Gamma\tau}$ to be real numbers. It is also useful to reparametrize, up to a constant, the $\hat{K}_{\tau}$ and $\hat{P_{\tau}}$ as

[TABLE]

where the $\lambda_{k\ell\tau}$ and $u_{\Gamma\tau}$ are a reparametrization of $\gamma_{k\ell\tau}$ and $\upsilon_{\Gamma\tau}$, and this form is used to derive the qubit parametrization (see Section III).

We first examine the $\mathcal{N}=1$ SPD. The G-type SPD for $\mathcal{N}=1$ is defined by $\hat{\rho}_{G1}=\hat{K}_{1}$, which is a non-interacting many-body density matrix. If $\hat{K}_{1}=|\Psi_{0}\rangle\langle\Psi_{0}|$ where $|\Psi_{0}\rangle$ is a single Slater determinant, then $\hat{\rho}_{G1}$ corresponds with the Hartree-Fock wave function. Conversely, the B-type SPD for $\mathcal{N}=1$ is defined as $\hat{\rho}_{B1}=\hat{P}_{1}$. Normally $\hat{\rho}_{B1}$ corresponds to a mixed state, but if $\hat{P}_{1}=|\Psi_{at}\rangle\langle\Psi_{at}|$ with $|\Psi_{at}\rangle$ defined as a direct product of atomic states from all sites, then $\hat{\rho}_{B1}$ corresponds to a pure state.

The G-type SPD with $\mathcal{N}=2$ is defined as

[TABLE]

If the center projector is taken as $\hat{K}_{2}=|\Psi_{0}\rangle\langle\Psi_{0}|$ with $|\Psi_{0}\rangle$ being a single Slater determinant, then $\hat{\rho}_{G2}=|\Psi_{G}\rangle\langle\Psi_{G}|$, where $|\Psi_{G}\rangle=\hat{P}_{1}|\Psi_{0}\rangle$ is the Gutzwiller wave function (GWF). The B-type SPD for $\mathcal{N}=2$ is defined as

[TABLE]

If the center projector is taken as $\hat{P}_{1}=|\Psi_{at}\rangle\langle\Psi_{at}|$ with $|\Psi_{at}\rangle$ defined as a direct product of atomic states from all sites, then $\hat{\rho}_{B2}=|\Psi_{B}\rangle\langle\Psi_{B}|$ , where $|\Psi_{B}\rangle=\exp\left(\sum_{k\ell}\gamma_{k\ell}\hat{n}_{k\ell}\right)|\Psi_{at}\rangle$ and we have previously referred to $|\Psi_{B}\rangle$ as the Bearyswil wave-function Cheng2021206402 ; Cheng2021195138 . However, it should be noted that $|\Psi_{B}\rangle$ is distinct from the original Bearyswil wave-function $|\Psi^{\prime}_{B}\rangle$ Baeriswyl19870 , which is defined as $|\Psi^{\prime}_{B}\rangle=\exp\left(\alpha\sum_{k\ell}\epsilon_{k\ell}\hat{n}_{k\ell}\right)|\Psi_{\infty}\rangle$, where $\alpha$ is a variational parameter and $|\Psi_{\infty}\rangle$ is the fully projected GWF. Though $|\Psi^{\prime}_{B}\rangle$ is technically a special case of a G-type $\mathcal{N}=3$ SPD, which can also be evaluated using the SCDA, it yields the same energy as $|\Psi_{B}\rangle$ for the insulating phase in $d=\infty$ for the multi-orbital Hubbard model (see Section VI.4.3 and Figure 3).

The G-type SPD for $\mathcal{N}=3$ combines the variational power of $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$, and is defined as

[TABLE]

If $\hat{K}_{2}$ corresponds to a single Slater determinant $|\Psi_{0}\rangle$, then $\hat{\rho}_{G3}=|\Psi_{GB}\rangle\langle\Psi_{GB}|$, where $|\Psi_{GB}\rangle=\exp\left(\sum_{k\ell}\gamma_{k\ell}\hat{n}_{k\ell}\right)\hat{P}_{1}|\Psi_{0}\rangle$ is a mild generalization of the original Gutzwiller-Baeriswyl wave function $|\Psi^{\prime}_{GB}\rangle$ introduced by Otsuka (see Eq. (2.3) in Ref. Otsuka19921645 ), defined as $|\Psi^{\prime}_{GB}\rangle=\exp\left(\alpha\sum_{k\ell}\epsilon_{k\ell}\hat{n}_{k\ell}\right)\hat{P_{1}}|\Psi_{0}\rangle$, where $\alpha$ is a variational parameter. Finally, the $\mathcal{N}=3$ B-type SPD is defined as $\hat{\rho}_{B3}=\hat{P}_{1}\hat{K}_{2}\hat{P}_{2}\hat{K}_{2}^{\dagger}\hat{P}_{1}^{\dagger}$. If $\hat{P}_{2}=|\Psi_{at}\rangle\langle\Psi_{at}|$, then $\hat{\rho}_{B3}=|\Psi_{BG}\rangle\langle\Psi_{BG}|$ where $|\Psi_{BG}\rangle=\hat{P}_{1}\exp\left(\sum_{k\ell}\gamma_{k\ell}\hat{n}_{k\ell}\right)|\Psi_{at}\rangle$ and we have previously referred to $|\Psi_{BG}\rangle$ as the Baeriswyl-Gutzwiller wave-function. However, it should be noted that $|\Psi_{BG}\rangle$ is distinct from the original Baeriswyl-Gutzwiller wave-function $|\Psi_{BG}^{\prime}\rangle$ Dzierzawa19951993 , which is defined as $|\Psi^{\prime}_{BG}\rangle=\hat{P}_{1}\exp\left(\alpha\sum_{k\ell}\epsilon_{k\ell}\hat{n}_{k\ell}\right)|\Psi_{\infty}\rangle$, where $\alpha$ is a variational parameter and $|\Psi_{\infty}\rangle$ is the fully projected GWF. The $|\Psi^{\prime}_{BG}\rangle$ is technically a special case of a G-type $\mathcal{N}=4$ SPD, which can also be evaluated using the SCDA.

It is useful to understand how the SPD recovers the exact solution with increasing $\mathcal{N}$ Cheng2021195138 . The Trotter-Suzuki decomposition Suzuki1976183 , given as

[TABLE]

can be considered as a special case of the SPD where the variational parameters are set by the parameters of the Hamiltonian and the temperature. Given that the large $\mathcal{N}$ Trotter-Suzuki decomposition recovers the exact finite temperature density matrix, the SPD in the large $\mathcal{N}$ limit possesses sufficient variational power to solve the Hamiltonian exactly. However, for a finite $\mathcal{N}$, the SPD leverages the full variational freedom while maintaining the integer time structure, leading to an observed exponential decrease of the energy error with increasing $\mathcal{N}$ Cheng2021206402 ; Cheng2022205129 . In contrast, the energy error for the Trotter-Suzuki decomposition diminishes polynomially with increasing $\mathcal{N}$ Suzuki1976183 . This key difference is the basis for the accuracy of the SPD for small $\mathcal{N}$.

Our initial efforts using VDAT have been exactly evaluating the SPD in $d=\infty$ via the SCDA, and this same approach can be used as a local approximation in finite dimensions. However, there is an existing body of literature which used variational quantum Monte-Carlo (VMC) to evaluate a simplified SPD in the Hubbard model. In the case of $\hat{\rho}_{G2}$, it is most convenient to execute VMC in real space, which has been used to evaluate the Gutzwiller wave function in the Hubbard model Yokoyama19871490 ; Yokoyama19873582 ; Yokoyama19882482 . For $\mathcal{N}>2$ in the one and two dimensional Hubbard model, a restricted form of the SPD has been evaluated using quantum Monte-Carlo to solve the single-orbital Hubbard model in two dimensions Otsuka19921645 ; Yamaji1998225 ; Yanagisawa19983867 ; Yanagisawa19993608 ; Koike200365 ; Yanagisawa2016114707 ; Yanagisawa2019054702 ; Yanagisawa20202040046 ; Yanagisawa202112 ; Yanagisawa2021127382 ; Sorella2023115133 ; Levy2024013237 , the $p$-$d$ model Yanagisawa202127004 ; Yanagisawa2001184509 , and selected molecules Chen20234484 .

II.2 Gauge Symmetry of the SPD

A important aspect of the SPD is the corresponding gauge symmetry Cheng2022205129 ; Cheng2023035127 , which will prove to be an important tool for reducing the redundancy of the parametrization of the SPD. Consider the general SPD $\hat{\rho}=\hat{K}_{1}\hat{P}_{1}\dots\hat{K}_{\mathcal{N}}\hat{P}_{\mathcal{N}}$, where $\hat{K}_{\tau}=\exp\left(\bm{\gamma}_{\tau}\cdot\hat{\bm{n}}\right)$ is a kinetic projector and $\boldsymbol{\gamma}\cdot\hat{\boldsymbol{n}}\equiv\sum_{\ell\ell^{\prime}}[\boldsymbol{\gamma}]_{\ell\ell^{\prime}}[\hat{\boldsymbol{n}}]_{\ell\ell^{\prime}}$ and $[\hat{\boldsymbol{n}}]_{\ell\ell^{\prime}}=\hat{a}_{\ell}^{\dagger}\hat{a}_{\ell^{\prime}}$, and $\hat{P}_{\tau}$ is a general interacting projector. The essence of the gauge symmetry lies in the fact that only the total product $\hat{\rho}$ determines physical observables, while there is considerable flexibility in the decomposition of $\hat{\rho}$ into different time steps and the separation between kinetic and interacting projectors. In order to analyze the gauge symmetry, we utilize the discrete action theory Cheng2021195138 , which is a formalism for studying integer time correlation functions of the SPD. The discrete action theory can be seen as a generalization of the many-body Green’s function formalism to the case of integer time, including an integer time Dyson equation given as

[TABLE]

where \text{\bm{g}}_{0} is the non-interacting integer time Green’s functions, $\bm{g}$ is the interacting integer time Green’s functions, and $\boldsymbol{S}=\exp(-\boldsymbol{\Sigma}^{T})$ is the exponential form of the integer time self-energy, where $\boldsymbol{\Sigma}$ is the integer time self-energy. The integer time Dyson equation recovers the usual Dyson when the SPD is chosen as the Trotter-Suzuki decomposition in the large $\mathcal{N}$ limit Cheng2021195138 . The \text{\bm{g}}_{0} and $\bm{g}$ are naturally defined in the compound space Cheng2021195138 as \text{\bm{g}}_{0}=\left\langle\textrm{\@text@baccent{\hat{\bm{n}}}}\right\rangle_{\textrm{\@text@baccent{\hat{\spdsymb}}}_{0}} and \text{\bm{g}}=\left\langle\textrm{\@text@baccent{\hat{\bm{n}}}}\right\rangle_{\textrm{\@text@baccent{\hat{\spdsymb}}}}, where [\textrm{\@text@baccent{\hat{\bm{n}}}}]_{\ell\tau,\ell^{\prime}\tau^{\prime}}=\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\dagger(\tau)}\textrm{\@text@baccent{\hat{a}}}_{\ell^{\prime}}^{(\tau^{\prime})}, the underbar indicates an operator in the compound space, $\ell=1,\dots,L$ enumerates the spin-orbitals for the entire system, and $\tau=1,\dots,\mathcal{N}$ labels the integer time step. The discrete Dyson equation is a matrix equation of dimension of $L\mathcal{N}\times L\mathcal{N}$, and it exactly relates the interacting and noninteracting integer time Green’s function via $\boldsymbol{S}$. Given that $\boldsymbol{S}$ is more convenient to work with, we refer to $\boldsymbol{S}$ as the integer time self-energy for brevity. For the convenience of discussing the gauge transformation of the SPD, it is useful to define several additional quantities

[TABLE]

where $[\boldsymbol{S}_{Q}]_{\tau\tau^{\prime}}=(-\delta_{\tau+1,\tau^{\prime}}+\delta_{\tau-\mathcal{N}+1,\tau^{\prime}})\boldsymbol{1}$ and $\left[\bm{S}_{K}\right]_{\tau\tau^{\prime}}=\delta_{\tau\tau^{\prime}}\exp\left(-\bm{\gamma}_{\tau}^{T}\right)$ are $L\times L$ matrices. Using these definitions, $\bm{g}_{0}$ and $\bm{g}$ can be defined succinctly as

[TABLE]

Generally, there are two types of gauge transformations of the SPD: intra-time-step and inter-time-step. The intra-time-step transformation occurs at the boundary of the kinetic and interacting projector for time step $\tau$ in the following way: $\hat{K}_{\tau}\rightarrow\hat{K}_{\tau}\exp\left(\bm{\mu}_{\tau}\cdot\hat{\bm{n}}\right)$ and $\hat{P}_{\tau}\rightarrow\exp\left(-\bm{\mu}_{\tau}\cdot\hat{\bm{n}}\right)\hat{P}_{\tau}$. Since the intra-time-step transformation does not alter the measurement at integer times, $\bm{g}$ is invariant to this transformation. However, given that $\hat{K}_{\tau}$ changes after the intra-time-step transformation, we have $\bm{S}^{\prime}_{K}=\bm{S}_{K}\bm{N}_{a}$, where $\left[\bm{N}_{a}\right]_{\tau\tau^{\prime}}=\delta_{\tau\tau^{\prime}}\exp\left(-\bm{\mu}_{\tau}^{T}\right)$, while $\bm{S}^{\prime}_{F}=\bm{S}_{F}$. Conversely, the inter-time-step transformation is defined at the boundary of the interacting projector at $\tau$ and the kinectic projector at $\tau+1$ in the following way: $\hat{P}_{\tau}\rightarrow\hat{P}_{\tau}\exp\left(-\bm{\mu}_{\tau}\cdot\hat{\bm{n}}\right)$ and $\hat{K}_{\tau+1}\rightarrow\exp\left(\bm{\mu}_{\tau}\cdot\hat{\bm{n}}\right)\hat{K}_{\tau+1}$, where $\bm{\mu}_{\mathcal{N}}\equiv 0$. We first discuss how $\bm{g}$ changes under this inter-time-step transformation, and it is useful to note the following identities (for details, see Eqns. (A22) and (A23) in Ref. Cheng2021195138 and Eqns. (A1)-(A4) in Ref. Cheng2022205129 )

[TABLE]

Using the above two equations, we can determine the transformation $\bm{g}^{\prime}=\bm{N}_{b}\bm{g}\bm{N}_{b}^{-1}$, where $\left[\bm{N}_{b}\right]_{\tau\tau^{\prime}}=\delta_{\tau\tau^{\prime}}\exp\left(-\bm{\mu}_{\tau}^{T}\right)$, and correspondingly $\bm{S}^{\prime}_{F}=\bm{N}_{b}\bm{S}_{F}\bm{N}_{b}^{-1}$. Using the transformation for the kinetic projector, we have $\bm{S}^{\prime}_{K}=\tilde{\bm{N}}_{b}\bm{S}_{K}$, where $\left[\tilde{\bm{N}}_{b}\right]_{\tau\tau^{\prime}}=\delta_{\tau\tau^{\prime}}\exp\left(-\bm{\mu}_{\tau-1}^{T}\right)$ and $\bm{\bm{\mu}}_{0}\equiv 0$, and $\tilde{\bm{N}}_{b}$ can be rewritten as $\tilde{\bm{N}}_{b}=\bm{S}_{Q}^{-1}\bm{N}_{b}\bm{S}_{Q}$. It should be noted that a general gauge transformation includes both intra-time-step and inter-time-step transformations, which can be parametrized by $\bm{N}_{a}$ and $\bm{N}_{b}$, and the above results can be synthesized as

[TABLE]

Moreover, the form of the integer time Dyson equation is invariant after the transformation, which requires

[TABLE]

Using Eq. (20), we have

[TABLE]

while using Eqns. (23), (19), (21), and (22), we have

[TABLE]

which can also be directly verified by assuming $\hat{P}_{\tau}$ is non-interacting. Correspondingly, the non-interacting integer time Green’s function changes as

[TABLE]

Notice that Eqns. (18)-(25) are exact relations. Within the SCDA Cheng2021195138 , $\bm{\mathcal{G}}$ will transform similarly to $\bm{g}_{0}$, where $\bm{N}_{a}$ and $\bm{N}_{b}$ are local. These gauge transformations provide a rigorous theoretical foundation for simplifying the task of achieving self-consistency within the SCDA. Moreover, choosing the appropriate gauge condition will also help to improve the numerical stability when minimizing over the variational parameters.

II.3 Review of the SCDA

The discrete action theory is a formalism designed to study integer time correlation functions of the SPD, and the self-consistent canonical discrete action theory (SCDA) is the integer time analogue of the dynamical mean-field theory Cheng2021206402 ; Cheng2021195138 . The SCDA exactly evaluates the SPD in infinite dimensions, and can be used as a local approximation for the SPD in finite dimensions. To appreciate the idea of the SCDA, we define the discrete action which encodes the integer time correlation of the SPD as

[TABLE]

where \textrm{\@text@baccent{\hat{\spdsymb}}}_{0}=\textrm{\@text@baccent{\hat{Q}}}\prod_{\tau=1}^{\mathcal{N}}\textrm{\@text@baccent{\hat{K}}}_{\tau}^{\left(\tau\right)} and \textrm{\@text@baccent{\hat{P}}}=\prod_{\tau=1}^{\mathcal{N}}\textrm{\@text@baccent{\hat{P}}}_{\tau}^{\left(\tau\right)} are the non-interacting and interacting parts of the discrete action ̱\hat{\varrho}, and the underbar denotes operators in the compound space.

The SCDA can be formulated as a scheme to evaluate the kinetic energy and local interacting energy through two self-consistent approximations of ̱\hat{\varrho} Cheng2022205129 , denoted as \textrm{\@text@baccent{\hat{\rho}}}_{K} and \textrm{\@text@baccent{\hat{\rho}}}_{loc}, respectively. Within the SCDA, the total energy is given as

[TABLE]

where \textrm{\@text@baccent{\hat{\rho}}}_{K} and \textrm{\@text@baccent{\hat{\rho}}}_{loc} are defined as

[TABLE]

where $\bm{\mathcal{G}}$ is the non-interacting integer time Green’s function for the local canonical discrete action \textrm{\@text@baccent{\hat{\rho}}}_{loc}. Both $\bm{S}$ and $\bm{\mathcal{G}}$ are $L\mathcal{N}\times L\mathcal{N}$ matrices, where $L$ is the number of spin-orbitals of the entire system. In the SCDA, the key assumption is that $\bm{S}$ and $\bm{\mathcal{G}}$ are block diagonal in real space clusters, given by $\left[\bm{S}\right]_{ij}=\bm{S}_{i}\delta_{ij}$ and $\left[\bm{\mathcal{G}}\right]_{ij}=\bm{\mathcal{G}}_{i}\delta_{ij}$ where $i$ and $j$ are real space cluster indices. To determine $\bm{S}_{i}$ and $\bm{\mathcal{G}}_{i}$, we have the following two sets of conditions for each cluster index $i$,

The integer time Dyson equation, given as

[TABLE]

provides a way to compute $\bm{S}_{i}$ by evaluating \textrm{\@text@baccent{\hat{\rho}}}_{loc;i}=\text{Tr}_{/i}\textrm{\@text@baccent{\hat{\rho}}}_{loc}, which can be expressed as

[TABLE] 2. 2.

The self-consistency condition for the local integer time Green’s function, given as

[TABLE]

where \boldsymbol{g}_{i}=\left\langle\textrm{\@text@baccent{\hat{\boldsymbol{n}}}}_{i}\right\rangle_{\textrm{\@text@baccent{\hat{\rho}}}_{loc}} and \boldsymbol{g}^{\prime}_{i}=\left\langle\textrm{\@text@baccent{\hat{\boldsymbol{n}}}}_{i}\right\rangle_{\textrm{\@text@baccent{\hat{\rho}}}_{K}}.

It is interesting to examine the sufficiency of Eq. (30) and (32) to determine both $\bm{\mathcal{G}}$ and $\bm{S}$ by counting the number of unknown entries and equations. We denote the number of sites in the lattice as $N_{site}$ and the number of spin orbitals on site $i$ as $N_{i}$. The number of unknowns in $\bm{\mathcal{G}}$ and $\bm{S}$ is $2\sum_{i=1}^{N_{site}}\left(N_{i}\mathcal{N}\right)^{2}$. For a given $i$, the two matrix equations given by Eq. (30) and (32) provide $2\left(N_{i}\mathcal{N}\right)^{2}$ entries. Therefore, we have a sufficient number of equations to determine $\bm{\mathcal{G}}$ and $\bm{S}$.

II.4 The tensor product representation for expectation values under \textrm{\@text@baccent{\hat{\rho}}}_{loc;i}

The key computational cost of implementing the SCDA is to compute local observables under \textrm{\@text@baccent{\hat{\rho}}}_{loc;i}, and a convenient approach was devised for $\mathcal{N}=3$ using a tensor product representation Cheng2023035127 . Here we generalize the tensor product representation to arbitrary $\mathcal{N}$. We assume that the interacting projector has the form $\hat{P}_{\tau;i}=\sum_{\Gamma}u_{\tau i\Gamma}\hat{X}_{i\Gamma}$ and $\boldsymbol{\mathcal{G}}_{i}$ is diagonal in the spin orbital index $\ell$, given as $\left[\bm{\mathcal{G}}_{i}\right]_{\ell\ell}=\delta_{\ell\ell^{\prime}}\mathcal{G}_{i,\ell}$, where $\mathcal{G}_{i,\ell}$ is an $\mathcal{N\times\mathcal{N}}$ matrix. For a local operator ̱\hat{O} at site $i$, the expectation value can be written as

[TABLE]

where the non-interacting discrete action for \textrm{\@text@baccent{\hat{\rho}}}_{loc;i,0} is given by

[TABLE]

where the vector $u_{\tau}=\left(u_{\tau i1},\dots,u_{\tau iN_{\Gamma}}\right)$ enumerates over the index $\Gamma$ with $N_{\Gamma}$ as the number of local diagonal Hubbard operators at site $i$, the $\mathcal{N}$ dimensional tensor \left(\textrm{\@text@baccent{\hat{O}}}\right)_{u} is defined as

[TABLE]

the direct product is defined as $\left[\otimes_{\tau=1}^{\mathcal{N}}u_{\tau}\right]_{\Gamma_{1},\dots,\Gamma_{\mathcal{N}}}=\prod_{\tau=1}^{\mathcal{N}}u_{\tau i\Gamma_{\tau}}$, and $\bm{A}\cdot\bm{B}$ is the contraction between two $\mathcal{N}$-dimensional tensors as

[TABLE]

We now consider operators that can be written as a product over the spin orbitals, given as \textrm{\@text@baccent{\hat{O}}}=\prod_{\ell}\textrm{\@text@baccent{\hat{O}}}_{\ell}. Assuming that $\mathcal{\bm{\mathcal{G}}}_{i}$ is diagonal in the spin orbital index $\ell$, we have

[TABLE]

where $\hat{X}{}_{i\Gamma;\ell}$ is defined in Eq. (4). Therefore, we can write \left(\textrm{\@text@baccent{\hat{O}}}\right)_{u} as a direct product

[TABLE]

where the direct product is defined as

[TABLE]

with \left(\textrm{\@text@baccent{\hat{O}}}_{\ell}\right)_{u;\ell} defined as

[TABLE]

which only depends $\left[\bm{\mathcal{G}}_{i}\right]_{\ell\ell}$. For some cases, there is no interacting projector on a given time step, and the corresponding dimension of \left(\textrm{\@text@baccent{\hat{O}}}_{\ell}\right)_{u;\ell} can be traced out. Moreover, Eq. (33) can be conveniently rewritten into a matrix product for all cases we study, which will be further discussed in Sections IV, V, and VI.

III Overview of the qubit energy form

The SCDA can exactly evaluate an SPD in $d=\infty$ with arbitrary $\mathcal{N}$ via a self-consistent solution, as outlined in Section II, and therefore the total energy can only be numerically evaluated for a given set of variational parameters in general. However, the recently developed gauge constrained SCDA algorithm Cheng2023035127 allows for an explicit evaluation of a G-type SPD with $\mathcal{N}\leq 3$, circumventing the need for self-consistent solution (see Section VI.2 for a review). In the present work, we offer a further refinement of the gauge constrained SCDA by transforming all variational parameters into physically intuitive variables, referred to as the qubit parametrization. Moreover, the qubit parametrization analytically resolves one constraint, reducing the number of variational parameters by one per spin orbital. Given the complexity of deriving the qubit parametrization, we collect all key results in this section, providing a self-contained presentation of all details needed to implement a practical calculation. A G-type and B-type SPD at $\mathcal{N}$ will be denoted as $\hat{\rho}_{G\mathcal{N}}$ and $\hat{\rho}_{B\mathcal{N}}$, respectively. For pedagogical purposes, we first present results for $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$ before finally considering $\hat{\rho}_{G3}$, illustrating that $\hat{\rho}_{G3}$ can be seen as a unification of $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$. Detailed derivations of all results for $\hat{\rho}_{G2}$, $\hat{\rho}_{B2}$, and $\hat{\rho}_{G3}$ are provided in Sections IV, V, and VI, respectively, while explicit applications are presented in Section VII. It should be emphasized that all qubit parametrization results are mathematically identical to corresponding VDAT results.

For the sake of generality, we consider the multiorbital Hubbard Hamiltonian given as

[TABLE]

where $H_{loc}$ is a polynomial function of the local density operators, given as

[TABLE]

where $c_{\ell_{1}\dots\ell_{n}}$ are parameters that define the n-body interactions. For a typical Hubbard model, only two-body interactions will be needed.

We begin by defining the local Hilbert space and the corresponding fermionic operators, and the corresponding map to qubit operators via the Jordan-Wigner transformation. For a given site $i$, the atomic configurations $\{|\Gamma\rangle\}$ form a complete basis for the local Hilbert space with dimension $2^{2N_{orb}}$, and we use $\hat{a}_{\ell}^{\dagger}$ and $\hat{a}_{\ell}$ to represent the creation and annihilation operators for spin-orbital $\ell$. Then we express the fermionic operators in terms of spin operators via the Jordan-Wigner transformation defined as

[TABLE]

where $\hat{\sigma}_{\ell}^{\pm}=\left(\hat{\sigma}_{\ell}^{x}\pm i\hat{\sigma}_{\ell}^{y}\right)/2$ and $\hat{\sigma}_{\ell}^{\mu}$ is defined as

[TABLE]

where $\mu\in\{x,y,z\}$ and $\hat{\sigma}^{\mu}$ is a standard $2\times 2$ Pauli matrix. It should be noted that we choose a convention for the Jordan-Wigner transformation such that $|0\rangle$ is associated with $|\uparrow\rangle$ and $|1\rangle$ is associated with $|\downarrow\rangle$. This local Jordan-Wigner transformation can be used to transform the usual Gutzwiller energy form into the qubit energy form, though spin operators will naturally emerge when using the tensor representation of local observables within the SCDA.

We proceed by evaluating the ground state energy under $\hat{\rho}_{G2}$, where $\hat{\rho}_{G2}$ is parametrized as

[TABLE]

where the variational parameters $\{u_{\Gamma}\}$ and $\{n_{k\ell,0}\}$ are real and $n_{k\ell,0}\in[0,1]$. The variational parameters $\{u_{\Gamma}\}$ can be reparametrized as a pure state of an effective $2N_{orb}$ qubit system (i.e. a $2N_{orb}$ spin $1/2$ system) characterized by a many-body density matrix $\bm{\rho}$. The total trial energy $E(\bm{\rho},\{n_{k\ell,0}\})$ for $\hat{\rho}_{G2}$ can then be written as

[TABLE]

with a constraint for each $\ell$

[TABLE]

Here we have taken the continuum limit of the discretized $n_{k\ell}$ and choose the convention $\int dk=1$. Equations (50)-(54) give an explicit functional form for the trial energy as a function of the variational parameters. The qubit energy form in Eqns. (50)-(54) is equivalent to a previous result obtained using the slave spin mean-field method De'medici2005205124 (see Appendix B for a detailed comparison). It should be emphasized that Eqns. (50)-(54) are simply a transformation of the SCDA algorithm, and therefore are completely equivalent to the usual Gutzwiller approximation.

We now proceed to study $\hat{\rho}_{B2}$, which is the dual of $\hat{\rho}_{G2}$, given as

[TABLE]

where the variational parameters $\{u_{\Gamma}\}$ and $\{\lambda_{k\ell}\}$ are real and $u_{\Gamma}\geq 0$. The qubit parametrization reparametrizes the $\{u_{\Gamma}\}$ and $\{\lambda_{k\ell}\}$ into $\bm{\rho}$ and $n_{k\ell}\in[0,1]$, where $\bm{\rho}$ is a diagonal, positive semi-definite $2^{2N_{orb}}\times 2^{2N_{orb}}$ matrix and $n_{k\ell}$ is the physical single particle density matrix, yielding the following total energy:

[TABLE]

with a constraint for each $\ell$

[TABLE]

It should be noted that the operators $\hat{n}_{eff,\ell}$ and $\hat{n}_{\ell}$ have the same expectation values under $\bm{\rho}$, which is evident from Eq. (57). This energy form provides a minimal description for the Mott insulating state in the multiorbital Hubbard model in $d=\infty$, much like the result of the Gutzwiller approximation for the metallic phase.

We now proceed to study $\hat{\rho}_{G3}$, which combines the variational capacity of both $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$. The key result of this paper is to recast the previously obtained explicit evaluation of the energy under $\hat{\rho}_{G3}$ into a physically intuitive form which can be viewed as a combination of the variational parameters of the qubit energy form for $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$. The $\hat{\rho}_{G3}$ is given as

[TABLE]

where the variational parameters $\{\lambda_{k\ell}\}$, $\{u_{\Gamma}\}$, and $\{n_{k\ell,0}\}$ are real and $n_{k\ell;0}\in[0,1]$. As demonstrated in Ref. Cheng2023035127 , when $\hat{K}\left(\{n_{k\ell,0}\}\right)$ corresponds to a Slater determinant (i.e. $n_{k\ell;0}=0$ or $n_{k\ell;0}=1$), the total energy can be explicitly evaluated, and therefore we follow this condition. Under this restricted form of $\{n_{k\ell,0}\}$, it is useful to introduce the concept of a reference Fermi surface, which delineates the boundary between $n_{k\ell;0}=0$ and $n_{k\ell;0}=1$ throughout the Brillouin zone. We can reparametrize $\{\lambda_{k\ell}\}$ and $\{u_{\Gamma}\}$ in terms of $\{n_{k\ell}\}$ and $\bm{\rho}$, where $\bm{\rho}$ is $2^{2N_{orb}}\times 2^{2N_{orb}}$ matrix corresponding to a pure state of a $2N_{orb}$ qubit system and $n_{k\ell}\in[0,1]$ is the physical momentum density distribution. The total trial energy is given as

[TABLE]

and there are two constraints for each $\ell$ given as

[TABLE]

where the symbol $<$ indicates that the integration is over the region where $n_{k\ell,0}=1$ for a given $\ell$, and the symbol $>$ indicates the region where $n_{k\ell,0}=0$. For a given $\ell$, the quantities $f_{\ell,0}$, $f_{\ell,x}$ and $f_{\ell,z}$ are nontrivial analytical functions of the following five variables

[TABLE]

where $\bm{\rho}$ and $n_{k\ell}$ are constrained such that

[TABLE]

which ensures that $f_{\ell,0}$, $f_{\ell,x}$ and $f_{\ell,z}$ are real numbers (see Eq. (73)). The explicit functional dependence is given by the following list of equations

[TABLE]

Several important points should be noted. First, the $\{n_{k\ell}\}$ only enter the local interaction energy through $\Delta_{\ell}$, $\mathcal{A}_{<\ell}$, $\mathcal{A}_{>\ell}$, and the constraint $\int dkn_{k\ell}=n_{\ell}$. Second, the $\{n_{k\ell,0}\}$ only enter the local interaction energy through the regions of integration (i.e. $>$ and $<$) and the constraint $\int dkn_{k\ell,0}=n_{\ell}$. Third, the operators $\hat{n}_{eff,\ell}$ and $\hat{n}_{\ell}$ have the same expectation values under $\bm{\rho}$, as in the case of the $\hat{\rho}_{B2}$. Finally, it should be clear that evaluating the $\hat{\rho}_{G3}$ has a similar computational cost as compared to $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$, where the largest computational cost is associated with evaluating the local interaction energy.

In the following, we briefly discuss how to numerically minimize the qubit energy form for $\hat{\rho}_{G3}$. For a multiorbital Hubbard model with $2N_{orb}$ spin-orbitals, the formal variational parameters are $\left\{n_{k\ell}\right\}$, $\left\{n_{k\ell;0}\right\}$, and $\bm{\rho}$ with $\ell=1,\dots,2N_{orb}$, and there are two local density constraints per spin orbital (see Eqn. (64)). We choose $\left\{n_{k\ell;0}\right\}$ corresponding to the momentum density distribution for the non-interacting Hamiltonian with local density $\left\{n_{\ell}\right\}$, and therefore $\left\{n_{k\ell;0}\right\}$ is determined from $\left\{n_{\ell}\right\}$. It is important to realize that the local interaction energy $\left\langle H_{loc}\left(\left\{\hat{n}_{eff,\ell}\right\}\right)\right\rangle_{\bm{\rho}}$ does not depend on the full details of $\{n_{k\ell}\}$, but instead only on $\left\{n_{\ell}\right\}$, $\left\{\Delta_{\ell}\right\},$ $\left\{\mathcal{A}_{<\ell}\right\}$, and $\left\{\mathcal{A}_{>\ell}\right\}$. Therefore, for a given spin orbital $\ell$, four Lagrange multipliers, $a_{<\ell}$, $b_{<\ell}$, $a_{>\ell}$, and $b_{>\ell}$, can be used to obtain the partially optimized $n_{k\ell}$ as Cheng2023035127

[TABLE]

where $X$ is either $<$ or $>$. Interestingly, the optimized $n_{k\ell}$ yields a non-trivial dependence on $\epsilon_{k\ell}$, in contrast to the Gutzwiller approximation (see Figure 1). There is a local density constraint between $\{n_{k\ell}\}$ and $\bm{\rho}$ given by Eqn. (64), and there are two strategies to enforce this. The first strategy is to start from $\{n_{k\ell}\}$, which are parametrized by $5\times 2N_{orb}$ variational parameters $\{n_{\ell}\}$, $\{a_{<\ell}\}$, $\{a_{>\ell}\}$, $\{b_{<\ell}\}$ and $\{b_{>\ell}\}$, though only $4\times 2N_{orb}$ are independent given Eqn. (64). Subsequently, the $\bm{\rho}$ can be parametrized by $2^{2N_{orb}}-2N_{orb}-1$ independent variational parameters Cheng2023035127 . Therefore, there are a total of $2^{2N_{orb}}+3\times 2N_{orb}-1$ independent variational parameters. This strategy motivates the construction of a one-body reduced density matrix functional, as presented in our companion paper companion . The second strategy is to start from $\bm{\rho}$, which can be parametrized by $2^{2N_{orb}}-1$ independent variational parameters, which determines $\{n_{\ell}\}$. Subsequently, the $\{\Delta_{\ell}\}$ , $\{b_{<\ell}\}$, and $\{b_{>\ell}\}$ are chosen as independent variational parameters, which determines $\{a_{X\ell}\}$, yielding $2^{2N_{orb}}+3\times 2N_{orb}-1$ independent variational parameters. The second strategy allows for a trivial implementation of the inequality constraint among $\Delta_{\ell}$, $n_{\ell}$, and $\xi_{\ell}$ (see Eq. (386)), for a given $\ell$, which is very important in practice. In either case, one is left with minimizing over $2^{2N_{orb}}+3\times 2N_{orb}-1$ variational parameters, which may be achieved using a variety of standard approaches.

Finally, we summarize the applications that are studied in Section VII.2. Consider the multiorbital Hubbard model with density-density interactions given as

[TABLE]

where $\hat{O}_{1}=\sum_{\alpha}\delta\hat{n}_{\alpha\uparrow}\delta\hat{n}_{\alpha\downarrow}$, $\hat{O}_{2}=\sum_{\alpha<\beta,\sigma}\delta\hat{n}_{\alpha\sigma}\delta\hat{n}_{\beta\bar{\sigma}}$, $\hat{O}_{3}=\sum_{\alpha<\beta,\sigma}\delta\hat{n}_{\alpha\sigma}\delta\hat{n}_{\beta\sigma}$, and $\delta\hat{n}_{\alpha\sigma}=\hat{n}_{\alpha\sigma}-\frac{1}{2}$, with the orbital indices $\alpha,\beta$ taking values of $1,\dots,N_{orb}$ and $\sigma\in\{\uparrow,\downarrow\}$. We consider the special case where $\epsilon_{k\ell}$ is independent of $\ell$ and has particle-hole symmetry. In this case, we show that Eq. (62) can be optimized over $\bm{\rho}$ and $\{n_{k\ell;0}\}$, with the assumption that the optimized $n_{k\ell;0}$ is given as $\theta(n_{k\ell}-n_{\ell}^{\star})$ where $\theta$ is the Heaviside function and $n_{\ell}^{\star}$ is chosen such that $\int dkn_{k\ell;0}=n_{\ell}$, yielding a trial energy purely as a functional of $\{n_{k\ell}\}$ as

[TABLE]

where $\widetilde{O}\left(\Delta\right)$ depends on $N_{orb}$ and $J/U$, and $n_{k\ell}$ is independent of $\ell$ and has particle-hole symmetry. Here we have dropped the spin-orbital index $\ell$ (e. g. $\Delta_{\ell}\rightarrow\Delta$) and define $A=\mathcal{A}_{<}+\mathcal{A}_{>}$. An important result is that $\widetilde{O}\left(\Delta\right)$ is non-analytic at $\Delta=\Delta_{c}$, which can be used to infer the nature of the Mott transition. For $N_{orb}>1$, the $\tilde{O}\left(\Delta\right)$ can be numerically represented by a one dimensional spline function companion , while $N_{orb}=1$ has an explicit expression which can be derived as

[TABLE]

The above equations finally deliver a closed form expression for a trial energy which qualitatively and quantitatively captures the Mott transition in the single band Hubbard model at half-filling in $d=\infty$, exactly reproducing the VDAT results for a G-type SPD at $\mathcal{N}=3$.

IV Derivation of the qubit energy form for $\hat{\rho}_{G2}$

In this section, we derive Eq. (50) using several different approaches. In Section IV.1, we review the GA using both a heuristic derivation and the central point expansion (CPE), and convert the standard energy form of the GA to the qubit energy form using the Jordan-Wigner transformation. In Section IV.2, we use the SCDA within the Gutzwiller gauge to evaluate $\hat{\rho}_{G2}$, and the tensor product representation is used to obtain the qubit energy form.

IV.1 The derivation of the qubit energy form from the GA

In this section, we provide an elementary derivation of the qubit energy form, which consists of two steps. In the first step, we derive the standard form of the GA, using both a heuristic argument and the CPE Cheng2020081105 , where the energy is parametrized by $\left\{n_{k\ell;0}\right\}$ and $\rho_{loc}$, where $\rho_{loc}$ is the local reduced density matrix of the ansatz, which is a diagonal matrix in the basis of $|\Gamma\rangle$. The quasi-particle weight is constructed from $\rho_{loc}$ and the fermionic creation and annihilation operators. In the second step, we convert the standard form of the GA to the qubit energy form which is parametrized by $\left\{n_{k\ell;0}\right\}$ and $\bm{\rho}$, where $\bm{\rho}$ is a pure state in a $2N_{orb}$-qubit space, and the quasi-particle weight is constructed from $\bm{\rho}$ and the Pauli spin operators.

IV.1.1 Derivation of the GA: a heuristic argument and the CPE

While the Gutzwiller approximation is well known, it is normally applied to the special case where $\hat{\rho}_{G2}$ is a pure state. Here we provide a heuristic derivation for a general $\hat{\rho}_{G2}$. Additionally, we use the CPE to derive the same result from a different perspective. We begin by presenting the heuristic derivation of the GA. First consider the expectation values for the diagonal Hubbard operator $\hat{X}_{i\Gamma}$ measured under $\hat{K}_{2}=\hat{K}\left(\left\{n_{k\ell;0}\right\}\right)$, given as

[TABLE]

where $n_{\ell}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2}}$ denotes the local density at a given site $i$ for the spin orbital $\ell$ under $\hat{K}_{2}$. The first assumption of the GA is that the atomic configuration distribution under $\hat{\rho}_{G2}$ can be approximated as

[TABLE]

which ignores off-site contributions in $\hat{K}_{2}$. Additionally, a constraint is applied to $\hat{P}_{1}$ such that the local density is invariant, given as

[TABLE]

which we refer to as the Gutzwiller constraint. The denominator in Eq. (102) ensures that the sum of all expectation values for the diagonal Hubbard operators are normalized. Given that the off-diagonal Hubbard operators have zero expectation value due to the restriction of the SPD used in this study, the expectation values for all diagonal Hubbard operators allow the evaluation of any local observable. The second assumption of the GA lies within the evaluation for the hopping term between two distinct sites. The idea is to break a hopping term like $\hat{a}_{i\ell}^{\dagger}\hat{a}_{j\ell}$ into two processes: first annihilating an electron at site $j$ and then creating an electron at site $i$, and the GA assumes the probabilities for the two steps are independent and only depend on the atomic distributions on site $i$ and $j$. The probability of creating or destroying an electron is renormalized by the interacting projector, and the GA assumes this is obtained by counting all relevant one-particle excitation processes as

[TABLE]

To simply this expression, we introduce $\rho_{loc}$ and $\rho_{loc;0}$, which are defined within the local Hilbert space as

[TABLE]

which are matrices of dimension $2^{2N_{orb}}\times 2^{2N_{orb}}$, yielding

[TABLE]

The Gutzwiller constraint, Eq. (103), can be rewritten as $\langle\hat{n}_{\ell}\rangle_{\rho_{loc}}=\langle\hat{n}_{\ell}\rangle_{\rho_{loc;0}}$. Given the second assumption of the GA, the single particle density matrix between two different sites $j$ and $j^{\prime}$ is renormalized as

[TABLE]

Combining Eq. (108) with the Gutzwiller constraint, the single particle density matrix for arbitrary $j$ and $j^{\prime}$ is given as

[TABLE]

The momentum density distribution $n_{k\ell}=\langle\hat{a}_{k\ell}^{\dagger}\hat{a}_{k\ell}\rangle_{\hat{\rho}_{G2}}$ can then be computed as

[TABLE]

where $N_{site}$ is the number of $k$-points, $\hat{a}_{k}^{\dagger}=\left(1/\sqrt{N_{site}}\right)\sum_{j}e^{ik\cdot j}\hat{a}_{j}^{\dagger}$, and $n_{k\ell;0}=\langle a_{k\ell}^{\dagger}a_{k\ell}\rangle_{\hat{K}_{2}}$. Therefore, the standard form for the total energy per site of the GA is given by

[TABLE]

with the density constraint

[TABLE]

The variational parameters satisfy $n_{kl;0}\in\left[0,1\right]$ and $\rho_{loc}$ is a diagonal positive semi-definite matrix.

We now provide an overview of how to derive Eq. (112) using the CPE (see Appendix C.1). We begin by making several observations about the GA. First, when $n_{k\ell;0}=n_{\ell}$, the $\hat{\rho}_{G2}$ describes a collection of atoms, and in this case Eq. (112) yields an exact evaluation in any dimension. Second, there is a linear relation between $n_{k\ell}$ and $n_{k\ell;0}$ given by

[TABLE]

Third, $\rho_{loc}$ is directly determined by $n_{\ell}$ and $\{u_{\Gamma}\}$, and is independent of $n_{k\ell;0}$. Recall that the $\hat{\rho}_{G2}$ is defined as $\hat{\rho}_{G2}=\hat{P}_{1}\hat{K}_{2}\hat{P}_{1}$, where $\hat{K}_{2}=\hat{K}\left(\left\{n_{k\ell;0}\right\}\right)$ and $\hat{P}_{1}=\hat{P}\left(\left\{u_{\Gamma}\right\}\right)$. We now introduce $\hat{K}_{2}^{\star}=\hat{K}\left(\left\{n_{k\ell;0}=n_{\ell}\right\}\right)$, which has the same local density matrix as $\hat{K}_{2}$, and $\rho_{G2}^{\star}=\hat{P}_{1}\hat{K}_{2}^{\star}\hat{P}_{1}$, which describes a collection of atoms. The $\{u_{\Gamma}\}$ can be reparametrized by the local reduced density matrix $\rho_{loc}^{\star}$ of $\rho_{G2}^{\star}$, where $\left[\rho_{loc}^{\star}\right]_{\Gamma\Gamma^{\prime}}=\delta_{\Gamma\Gamma^{\prime}}\langle\hat{X}_{i\Gamma}\rangle_{\rho_{G2}^{\star}}$ and is constrained by $n_{\ell}=\langle\hat{a}_{\ell}^{\dagger}\hat{a}_{\ell}\rangle_{\rho_{loc}^{\star}}.$ Therefore, $n_{k\ell}$ and $\rho_{loc}$ are functionals of $\left\{n_{k\ell;0}\right\}$ and $\rho_{loc}^{\star}$. The CPE amounts to the expansion of observables in terms of $\{n_{k\ell;0}\}$ about $\left\{n_{k\ell;0}=n_{\ell}\right\}$, and up to the first order, one recovers Eq. (116) and $\rho_{loc}=\rho_{loc}^{\star}$, proving that the first-order CPE recovers the GA.

IV.1.2 Converting the GA energy into the qubit form using the Jordan-Wigner

transformation

Here we discuss how to convert the standard form of the GA energy into the qubit form, which is mathematically equivalent. The qubit form provides a unified view of the energy evaluated using $\hat{\rho}_{G2}$, $\hat{\rho}_{B2}$, and $\hat{\rho}_{G3}$. The qubit parametrization consists of two steps. First, we introduce a purified many-body density matrix $\bm{\rho}$ from $\rho_{loc}$. Second, we perform the Jordan-Wigner transformation, which converts $\mathcal{R}_{\ell}$ (Eq. (114)) from an expression involving fermionic operators $\hat{a}_{\ell}^{\dagger}$ and $\hat{a}_{\ell}$ into an expression involving spin operators $\hat{\sigma}_{\ell}^{x}$. We begin by defining $\bm{\rho}=|\Psi\rangle\langle\Psi|$, where $|\Psi\rangle=\sum_{\Gamma}\sqrt{\left[\rho_{loc}\right]_{\Gamma\Gamma}}|\Gamma\rangle$, which yields

[TABLE]

Using the Jordan-Wigner transformation (see Eqns. (43) and (44)), we obtain

[TABLE]

Therefore, we can rewrite the numerator of $\mathcal{R}_{\ell}$ (see Eq. (107)) as

[TABLE]

where we used the fact that $\bm{\rho}$ is a real symmetric matrix and the following relation

[TABLE]

Moreover, the denominator of $\mathcal{R}_{\ell}$ (see Eq. (107)) is given as

[TABLE]

Therefore, we have demonstrated that $\mathcal{R}_{\ell}=\xi_{\ell}/\xi_{\ell;0}$, and the local interaction can be written as

[TABLE]

Using the expression for $\mathcal{R}_{\ell}$ and the local energy, we arrive at the energy expressions given in Eqns. (50)-(53). This proof demonstrates that the usual Gutzwiller approximation can be straightforwardly transformed into the qubit energy form, which is equivalent to the result of the slave spin mean-field theory De'medici2005205124 (see Appendix B for additional details).

IV.2 The derivation of the qubit energy form using the SCDA

In this section, we use the gauge constrained SCDA to derive the qubit energy form. The derivation consists of two steps. First, in Section IV.2.1, the Gutzwiller gauge is used to rederive the standard form of the GA energy. Second, in Section IV.2.2, we derive the qubit energy form using the tensor representation, which does not rely on the Jordan-Wigner transformation. Additionally, in Section IV.2.3, we discuss how the quantities $\boldsymbol{\mathcal{G}}$ , $\boldsymbol{S}$, and $\boldsymbol{g}$ change under a general gauge transformation.

IV.2.1 The SCDA within the Gutzwiller

gauge

In this section, we demonstrate how the gauge of the SPD can be used to automatically satisfy the SCDA self-consistency condition Cheng2021195138 , and we utilize some notation from the CPE. Interestingly, the Gutzwiller constraint of the GA can be used to define an appropriate gauge for the SPD, which we refer to as the Gutzwiller gauge. Starting from $\hat{\rho}_{G2}=\hat{P}_{1}\hat{K}_{2}\hat{P}_{1}$ with an arbitrary gauge, a gauge transformation $\hat{P}_{1}\rightarrow\hat{P}_{1}\hat{N}^{-1}$ and $\hat{K}_{2}\rightarrow\hat{N}\hat{K}_{2}\hat{N}$ can always be performed to ensure that $\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{\rho}_{G2}}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2}}=n_{\ell}$, where $\hat{N}=\exp\left(\sum_{i\ell}\mu_{\ell}\hat{n}_{i\ell}\right)$ is chosen to satisfy the Gutzwiller constraint. Within this Gutzwiller gauge, we can choose $\left[\bm{\mathcal{G}}_{i}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}\mathcal{G}_{\ell}$, where the component for the spin-orbital $\ell$ is given by

[TABLE]

We now prove that this $\mathcal{G}_{\ell}$ ensures that the SCDA self-consistency condition is automatically satisfied. The first step of the proof relies on the fact that under the Gutzwiller gauge \textrm{\@text@baccent{\hat{\rho}}}_{loc} within the SCDA is the discrete action of the central point of the SPD $\rho_{G2}^{\star}=\hat{P}_{1}\hat{K}_{2}^{\star}\hat{P}_{1}$, which can be shown as follows. The local reduced density matrix of $\rho_{G2}^{\star}$ at site $i$ is $\hat{\rho}_{G2;i}^{\star}=\hat{P}_{1;i}\hat{K}_{2;i}^{\star}\hat{P}_{1;i}$, which yields a discrete action \textrm{\@text@baccent{\hat{\rho}}}_{G2;i}^{\star}=\textrm{\@text@baccent{\hat{\rho}}}_{G2;i;0}^{\star}\textrm{\@text@baccent{\hat{P}}}_{1;i}^{\left(1\right)}\textrm{\@text@baccent{\hat{P}}}_{1;i}^{\left(2\right)}, where \textrm{\@text@baccent{\hat{\rho}}}_{G2;i;0}^{\star}=\textrm{\@text@baccent{\hat{Q}}}\textrm{\@text@baccent{\hat{K}}}_{2;i}^{\star}. Given that $\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2;i}^{\star}}=n_{\ell}$, the integer time Green’s function for \textrm{\@text@baccent{\hat{\rho}}}_{G2;i;0}^{\star} is given by Eq. (127), and correspondingly \textrm{\@text@baccent{\hat{\rho}}}_{G2;i}^{\star} is equivalent to the local discrete action \textrm{\@text@baccent{\hat{\rho}}}_{loc;i} for the SCDA, proving \textrm{\@text@baccent{\hat{\rho}}}_{G2}^{\star}=\textrm{\@text@baccent{\hat{\rho}}}_{loc}. Local observables within the SCDA can then be evaluated under the central point of the SPD as

[TABLE]

consistent with the relation derived from the CPE (see Eq. (595)). Similarly, the local interacting integer time Green’s function $\left[\bm{g}_{i}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}g_{\ell}$ can be evaluated as

[TABLE]

where

[TABLE]

and

[TABLE]

and the Gutzwiller gauge ensures that the diagonal part of $g_{\ell}$ is $n_{\ell}$. To connect $\mathcal{R}_{\ell,12}$ and $\mathcal{R}_{\ell,21}$ with $\mathcal{R}_{\ell}$, we use the following relations

[TABLE]

which yields $\mathcal{R}_{\ell,12}=\mathcal{R}_{\ell,21}=\mathcal{R}_{\ell}$. The integer time self-energy can then be computed as $\left[\bm{S}_{i}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}S_{\ell}$, using $S_{\ell}=\left(\mathcal{G}_{\ell}^{-1}-1\right)^{-1}\left(g_{\ell}^{-1}-1\right)$, which yields

[TABLE]

Therefore, the interacting integer time Green’s function for a given $k$-point is given by $\left[\bm{g}_{k}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}g_{k\ell}$, where $g_{k\ell}=\left(1+\left(g_{k\ell;0}^{-1}-1\right)S_{\ell}\right)^{-1}$, which yields

[TABLE]

and the momentum density distribution $n_{k\ell}=\left[g_{k\ell}\right]_{22}$ is in agreement with Eq. (111). Finally, the SCDA self-consistency can be verified as $\frac{1}{N_{site}}\sum_{k}g_{k\ell}=g_{\ell}.$

IV.2.2 Derivation of the qubit energy form

Here we show how to derive the qubit energy form, given in Eqns. (50)-(54), using the tensor product representation. We begin by evaluating the relevant observables under \textrm{\@text@baccent{\hat{\rho}}}_{loc}, where the tensor product representation in Eq. (33) simplifies to

[TABLE]

given that $u^{T}=u_{1}=u_{2}$. The relevant components needed to evaluate the local integer time Green’s function in the Gutzwiller gauge are given by

[TABLE]

A linear transformation $u=Vw$ can be introduced such that

[TABLE]

where

[TABLE]

We choose $V_{\ell}$ such that the identity operator in the $w$-representation is an identity matrix, which is given by

[TABLE]

and correspondingly, the components in the $w$-representation are given by

[TABLE]

To connect with the qubit energy form, we define the many-body density matrix corresponding to a pure state of the qubit system as

[TABLE]

where the renormalization factor $\mathcal{R}_{\ell}$ can now be rewritten as

[TABLE]

using Eqns. (167), (170), and (129). It should be emphasized that $\hat{\sigma}_{\ell}^{x}$ is a $2^{2N_{orb}}\times 2^{2N_{orb}}$ matrix, defined in Eq. (46). Additionally, the local expectation value of the interaction energy can be rewritten as

[TABLE]

using Eq. (173) and (161), where $\hat{n}_{\ell}$ is defined in Eq. (45).

IV.2.3 SCDA under a general gauge transformation

Here we discuss the SCDA within an arbitrary gauge, where $\hat{\rho}_{G2}=\hat{P}_{1}\hat{K}_{2}\hat{P}_{1}^{\dagger}$, as this will be important to understanding the gauge constrained algorithm for $\mathcal{N}=3$. A general gauge transformation which maintains a G-type form is defined as $\hat{P}_{1}\rightarrow\hat{P}_{1}\hat{N}^{-1}$, $\hat{K}_{2}\rightarrow\hat{N}\hat{K}_{2}\hat{N}^{\dagger}$ , $\hat{P}_{1}^{\dagger}\rightarrow\left(\hat{N}^{\dagger}\right)^{-1}\hat{P}_{1}^{\dagger}$, where $\hat{N}=\exp\left(\bm{\mu}\cdot\hat{\bm{n}}\right)$. Using the results of Section II.2, this gauge transformation can be decomposed in terms of an intra-time-step transformation given by $\bm{N}_{a}=\text{diag}\left(\boldsymbol{1},\exp\left(-\bm{\mu}^{*}\right)\right)$, and an inter-time-step transformation given by $\bm{N}_{b}=\text{diag}\left(\exp\left(-\bm{\mu}^{T}\right),\boldsymbol{1}\right)$. The corresponding transformations for the integer time Green’s functions can be found in Section II.2. In the following, we focus on the transformation with the form $\left[\bm{\mu}\right]_{i\ell,i^{\prime}\ell}=\delta_{ii^{\prime}}\delta_{\ell\ell^{\prime}}\mu_{\ell}$ where $\mu_{\ell}$ is a real number. Using Eq. (18), the transformation of $\left[\bm{g}_{i}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}g_{\ell}$ with the component $g_{\ell}$ is given as

[TABLE]

Using Eq. (24), the transformation for $\bm{S}_{i}$ is given by

[TABLE]

Finally, the transformation for $\bm{\mathcal{G}}_{i}$ is given by

[TABLE]

We proceed by exploring a particular choice of gauge for the SPD, referred to as the anti-symmetric gauge, which will motivate the gauge choice for the case of $\mathcal{N}=3$. The essence of the anti-symmetric gauge is to choose $\mu_{\ell}$ such that $\left[g^{\prime}_{\ell}\right]_{12}=-\left[g^{\prime}_{\ell}\right]_{21}$, which can be accomplished as

[TABLE]

Under the anti-symmetric gauge, we have

[TABLE]

The anti-symmetric gauge can also automatically satisfy the SCDA self-consistency condition, and be used to derive the qubit energy form.

V Derivation of the qubit energy form for $\hat{\rho}_{B2}$

In this section, we derive Eq. (56) using several different approaches. In Section V.1, we present two derivations: a heuristic approach and the central point expansion (CPE). In Section V.2, we use the SCDA to evaluate $\hat{\rho}_{B2}$, and the tensor product representation is used to obtain the qubit energy form.

V.1 Derivation of the qubit energy form: a heuristic

approach and the CPE

To begin, it is useful to rewrite the local interaction Hamiltonian given in Eq. (42) as

[TABLE]

where $\delta\hat{D}_{iI}=\prod_{\ell\in I}\delta\hat{n}_{i\ell}$, the density fluctuation is defined as $\delta\hat{n}_{i\ell}=\hat{n}_{i\ell}-n_{\ell}$, the index $I$ enumerates all possible subsets of the local spin orbitals, with $\delta\hat{D}_{iI}=1$ when $I=\left\{\right\}$, and the parameters $E_{I}$ reparametrize the coefficients in Eq. (42). Using the alternative form of the local interaction Hamiltonian given in Eq. (195), the corresponding qubit form of the trial energy is given as

[TABLE]

where the variational parameters $n_{k\ell}\in\left[0,1\right]$, the $\bm{\rho}$ is a many-body density matrix for a $2N_{orb}$ qubit system which is diagonal in the Pauli-Z basis, $\delta\hat{D}_{I}=\prod_{\ell\in I}(\hat{n}_{\ell}-n_{\ell})$, and $n_{\ell}=\int dkn_{k\ell}$. For a given $\ell$, there is constraint given by

[TABLE]

Below, we present two approaches for deriving Eq. (196).

A heuristic approach for deriving Eq. (196) is via the formal duality between $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$. First, in $\hat{\rho}_{G2}$, the center projector $\hat{K}_{2}$ is constrained to have the same local density as $\hat{\rho}_{G2}$, i.e., $\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2}}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{\rho}_{G2}}$. Similarly, the center projector $\hat{P}_{1}$ in $\hat{\rho}_{B2}$ is constrained to have the same local density as $\hat{\rho}_{B2}$, i.e., $\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{P}_{1}}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{\rho}_{B2}}$. Second, in $\hat{\rho}_{G2}$ (see Eq. (50)), the momentum density fluctuation $\delta n_{k\ell}=n_{k\ell}-n_{\ell}$ is renormalized from the bare momentum density fluctuation $\delta n_{k\ell;0}=n_{k\ell;0}-n_{\ell}$ with a factor $\mathcal{Z}=\xi_{\ell}^{2}/\xi_{\ell;0}^{2}$, which depends on $\bm{\rho}$. To the contrary, $\langle\delta\hat{D}_{iI}\rangle_{\hat{\rho}_{B2}}$ is renormalized from the reference value $\langle\delta\hat{D}_{iI}\rangle_{\hat{P}_{1}}$ with a factor $\prod_{\ell\in I}\mathcal{F_{\ell}}$, where $\mathcal{F}_{\ell}$ depends on $n_{k\ell}$. The expression for $\mathcal{F}_{\ell}$ can be determined using a counting scheme similar to the one used within the GA. We begin by transforming $\delta\hat{D}_{iI}$ into momentum space, which requires a summation over terms consisting of $N_{I}$ creation and $N_{I}$ annihilation operators, where $N_{I}$ is the number of spin orbitals in set $I$. Each creation or annihilation process for $k\ell$ will be scaled by

[TABLE]

where $\Gamma\in\{1,2\}$ enumerates the empty and occupied states for a given $k\ell$ and

[TABLE]

Using the infinite dimensional approximation where momentum conservation can be neglected, we obtain $\langle\delta\hat{D}_{iI}\rangle_{\hat{\rho}_{B2}}=\left(\prod_{\ell\in I}\mathcal{F}_{\ell}\right)\langle\delta\hat{D}_{iI}\rangle_{\hat{P}_{1}}$. By identifying $\bm{\rho}$ as the local reduced density matrix of $\hat{P}_{1}$, we obtain the form given in Eq. (197).

A more rigorous way to derive the factor in Eq. (197) is to use the CPE. We outline the derivation here, and the details can be found in Appendix C.2. The CPE for $\hat{\rho}_{B2}$ is a dual version to the CPE of $\hat{\rho}_{G2}$, where the latter is described in Section IV.1.1. We start by defining the central point as $\hat{\rho}_{B2}^{\star}=\hat{K}_{1}\hat{P}_{1}^{\star}\hat{K}_{1}$, where $\hat{P}_{1}^{\star}$ is a non-interacting projector chosen such that $n_{\ell}\equiv\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{P}_{1}}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{P}_{1}^{\star}}$. The kinetic projector $\hat{K}_{1}$ can be parametrized using $n_{k\ell}^{\star}=\langle\hat{n}_{k\ell}\rangle_{\hat{\rho}_{B2}^{\star}}$ with the constraint$\int dkn_{k\ell}^{\star}=n_{\ell},$ while $\hat{P}_{1}$ can be parametrized using $\delta D_{iI;0}=\langle\delta\hat{D}_{iI}\rangle_{\hat{P}_{1}}$ and $\{n_{\ell}\}$. Therefore, the observables under $\hat{\rho}_{B2}$ are functionals of $n_{k\ell}^{\star}$ and $\left\{\delta D_{iI;0}\right\}$. Performing a first order expansion in terms of $\left\{\delta D_{iI;0}\right\}$ about $\left\{\delta D_{iI;0}=0\right\}$, we obtain $n_{k\ell}=n_{k\ell}^{\star}$ and

[TABLE]

Finally, we explain how to obtain Eqns. (56)-(58) from Eq. (202). We begin by defining the effective density operator $\hat{n}_{eff,\ell}=n_{\ell}+\mathcal{F}_{\ell}\delta\hat{n}_{\ell}$ where $\delta\hat{n}_{\ell}=\hat{n}_{\ell}-n_{\ell}$, the corresponding fluctuation form $\delta\hat{n}_{eff,\ell}\equiv\hat{n}_{eff,\ell}-\hat{n}_{\ell}=\mathcal{F}_{\ell}\delta\hat{n}_{\ell}$, and a diagonal many-body density matrix $\bm{\rho}$ for a $2N_{orb}$ qubit system $\left[\bm{\rho}\right]_{\Gamma\Gamma^{\prime}}=\delta_{\Gamma\Gamma^{\prime}}\langle\hat{X}_{i\Gamma}\rangle_{\hat{P}_{1}}$. Using $\delta D_{iI,0}=\left\langle\prod_{i\in I}\delta\hat{n}_{\ell}\right\rangle_{\bm{\rho}}$ and Eq. (202), we have $\delta D_{iI}=\left\langle\prod_{i\in I}\delta\hat{n}_{eff,\ell}\right\rangle_{\bm{\rho}}$. Furthermore, $n_{\ell}=\left\langle\hat{n}_{\ell}\right\rangle_{\bm{\rho}}$ implies that $\left\langle\delta\hat{n}_{\ell}\right\rangle_{\bm{\rho}}=\left\langle\delta\hat{n}_{eff,\ell}\right\rangle_{\bm{\rho}}=0$. Therefore, the expectation value of $H_{loc}\left(\left\{\hat{n}_{i\ell}\right\}\right)$ is given by $\left\langle H_{loc}\left(\left\{\hat{n}_{eff,\ell}\right\}\right)\right\rangle_{\bm{\rho}}$, providing the connection to Eqns. (56)-(58).

V.2 The derivation of the qubit energy form via the SCDA

In this section, we use the gauge constrained SCDA to derive the qubit energy form. The derivation consists of two steps. First, in Section V.2.1, we use the Gutzwiller gauge to automatically satisfy the SCDA self-consistency condition. Second, in Section V.2.2, we derive the qubit energy form using the tensor representation. Additionally, in Section V.2.3, we discuss how the quantities $\boldsymbol{\mathcal{G}}$ , $\boldsymbol{S}$, and $\boldsymbol{g}$ change under a general gauge transformation.

V.2.1 SCDA within the Gutzwiller gauge

In this section, we use the gauge constrained SCDA to evaluate $\hat{\rho}_{B2}=\hat{K}_{1}\hat{P}_{1}\hat{K}_{1}$, where $\hat{P}_{1}=\exp\left(\sum_{i\Gamma}\upsilon_{\Gamma}\hat{X}_{i\Gamma}\right)$ and $\hat{K}_{1}=\exp\left(\sum_{k\ell}\gamma_{k\ell}\hat{n}_{k\ell}\right)$, recovering the results from Section V.1. The key observation is that for $\hat{\rho}_{B2}$, the interacting projectors only act on the first integer time step. Consider a gauge transformation for $\hat{\rho}_{B2}$ given as $\hat{K}_{1}\rightarrow\hat{K}_{1}\hat{N}$ and $\hat{P_{1}}\rightarrow\hat{N}^{-1}\hat{P}_{1}\hat{N}^{-1}$ where $\hat{N}=\exp\left(\sum_{i\ell}\mu_{\ell}\hat{n}_{i\ell}\right)$. We choose the gauge transformation such that $\bm{\mathcal{G}}_{i}$ is given as $\left[\bm{\mathcal{G}}_{i}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}\mathcal{G}_{\ell}$, where

[TABLE]

and $n_{\ell}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{\rho}_{B2;0}}$ is the local density for the non-interacting SPD $\hat{\rho}_{B2;0}=\hat{K}_{1}\hat{K_{1}}$. Under this gauge choice, $\hat{P}_{1}$ will ensure that $\left[\bm{g}_{i}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}g_{\ell}$ with $\left[g_{\ell}\right]_{11}=\left[\mathcal{G}_{\ell}\right]_{11}=n_{\ell}$. The integer time self-energy \left[\bm{S}_{i}\right]$${}_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}S_{\ell} can then be determined as

[TABLE]

The remaining entries of Eq. (203) should be determined from the self-consistency condition in Eq. (32). Given that $S_{\ell}$ is the identity matrix, the non-interacting and interacting integer time Green’s functions are given by $\left[\bm{g}_{k;0}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}g_{k\ell;0}$ and $\left[\bm{g}_{k}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}g_{k\ell}$ where

[TABLE]

and $n_{k\ell}$ is the momentum density distribution. We can then verify that the self-consistency condition $\bm{g}_{i}=\left(1/N_{\textrm{site}}\right)\sum_{k}\bm{g}_{k}$ is fulfilled, given that

[TABLE]

where $n_{\ell}=\left(1/N_{site}\right)\sum_{k}n_{k\ell}$ and the $A_{\ell}$ is defined as

[TABLE]

We analogously refer to this gauge as the Gutzwiller gauge given that $g_{\ell}=\mathcal{G}_{\ell}$, and therefore the diagonal elements of the two matrices are the same.

V.2.2 Derivation of the

qubit energy form

In this section, we use the tensor product representation to derive the qubit energy form, paralleling the procedure in Section IV.2.2. Given that we only have interacting projectors at the first integer time step, operators in the $u$-representation \left(\textrm{\@text@baccent{\hat{O}}}\right)_{u} can be reduced from a matrix to a vector, given as

[TABLE]

and the expectation value of ̱\hat{O} at site $i$ is given as

[TABLE]

where $u=u_{1}$ and the dot denotes the normal dot product between two vectors. Given the direct product structure of \left(\textrm{\@text@baccent{\hat{O}}}\right)_{u}, we only need to compute the component for a given spin-orbital. Here, we only list the relevant matrix elements needed to derive the $w$-representation, given as

[TABLE]

The $u$-representation and $w$-representation are related by $u=Vw$ such that

[TABLE]

and therefore

[TABLE]

When $V$ has a direct product form $V=V_{1}\otimes\dots\otimes V_{2N_{orb}},$ then

[TABLE]

and $V_{\ell}$ is chosen to obtain \left(\textrm{\@text@baccent{\hat{1}}}\right)_{w;\ell} as a vector of ones. The matrix elements are given as

[TABLE]

From the preceding equations, we have

[TABLE]

where

[TABLE]

and therefore we have

[TABLE]

We now proceed to reinterpret the local energy in terms of the qubit representation, identifying \left(\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\dagger\left(1\right)}\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\left(1\right)}\right)_{w}=\textrm{diag}(\hat{n}_{\ell}) and \left(\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\dagger\left(2\right)}\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\left(2\right)}\right)_{w}=\textrm{diag}(\hat{n}_{eff,\ell}), Eq. (228) becomes

[TABLE]

Furthermore, we define a diagonal many-body density matrix $\bm{\rho}$ with $\textrm{diag}(\bm{\rho})=w$, resulting in

[TABLE]

Therefore, the qubit energy form has been recovered.

V.2.3 SCDA under a general gauge transformation

Here we discuss the SCDA within an arbitrary gauge, analogous to Section IV.2.3. Given $\hat{\rho}_{B2}=\hat{K}_{1}\hat{P}_{1}\hat{K}_{1}^{\dagger}$, a general gauge transformation is given as: $\hat{K}_{1}\rightarrow\hat{K}_{1}\hat{N}$, $\hat{P}_{1}\rightarrow\hat{N}^{-1}\hat{P}_{1}\left(\hat{N}^{\dagger}\right)^{-1}$, $\hat{K}_{1}^{\dagger}\rightarrow\hat{N}^{\dagger}\hat{K}_{1}^{\dagger}$, where $\hat{N}=\exp\left(\bm{\mu}\cdot\hat{\bm{n}}\right)$. Therefore, we have $\bm{N}{}_{a}=\text{diag}\left(\exp\left(-\bm{\mu}^{T}\right),\boldsymbol{1}\right)$ and $\bm{N}_{b}=\text{diag}\left(\exp\left(-\bm{\mu}^{*}\right),\boldsymbol{1}\right)$. Assuming $\left[\bm{\mu}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}\mu_{\ell}$, where $\mu_{\ell}$ is real, we can obtain the transformation for $g_{\ell}$, $S_{\ell}$, and $\mathcal{G}{}_{\ell}$ as

[TABLE]

We can define an anti-symmetric gauge, similar to the $\hat{\rho}_{G2}$ case, where we choose $\left[\mathcal{G}_{\ell}^{\prime}\right]_{11}=1/2$, leading to $\mu_{\ell}=\frac{1}{2}\ln\left(\frac{1-n_{\ell}}{n_{\ell}}\right)$, yielding

[TABLE]

The anti-symmetric gauge can also automatically satisfy the SCDA self-consistency condition, and be used to derive the qubit energy form.

VI Derivation of the qubit energy form for

$\hat{\rho}_{G3}$

In this section, we derive Eq. (62) using the gauge constrained SCDA. In section VI.1, we provide a high level comparison of the original gauge constrained algorithm and the qubit parametrization. In section VI.2, we provide a review of the original gauge constrained algorithm, which is necessary to understand the qubit parametrization in this work. In section VI.3, we propose the qubit parametrization which yields the qubit energy form. In section VI.4, we examine the case of half-filling, and explore how the qubit energy form for $\hat{\rho}_{G3}$ can recover the cases of $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$.

VI.1 Comparing the original gauge

constrained trial energy to the qubit energy form

In this section, we outline how the qubit parametrization improves the original gauge constrained algorithm Cheng2023035127 . In the original gauge constrained algorithm, the trial energy under $\hat{\rho}_{G3}$ is given as

[TABLE]

which is a function of $\{n_{k\ell}\}$, $\left\{\mathcal{G}_{12,\ell}\right\}$, $\{u_{\Gamma}\}$, and $\{n_{k\ell;0}\}$, given that $\left\{\mathcal{A}_{<\ell}\right\}$ and $\left\{\mathcal{A}_{>\ell}\right\}$ are functions of $\{n_{k\ell}\}$ and $\{n_{k\ell;0}\}$. It should be noted that $n_{k\ell;0}\in\{0,1\}$ and $n_{k\ell}\in[0,1]$, and there are three constraints for a given $\ell$

[TABLE]

where the functions $n_{\ell}\left(\left\{\mathcal{G}_{12;\ell}\right\},\{u_{\Gamma}\}\right)$ and $\Delta_{\ell}\left(\left\{\mathcal{G}_{12;\ell}\right\},\{u_{\Gamma}\}\right)$ are explicitly defined in Ref. Cheng2023035127 . While this parametrization allows for an explicit evaluation of the total energy, there are several shortcomings of this parametrization. First, the function $n_{\ell}\left(\left\{\mathcal{G}_{12;\ell}\right\},\{u_{\Gamma}\}\right)$ is highly non-trivial and thus the minimization under a fixed density is cumbersome. Previously, this problem was addressed by introducing a linear transformation over $\{u_{\Gamma}\}$, known as the $w$ representation. Second, the function $\Delta_{\ell}\left(\left\{\mathcal{G}_{12;\ell}\right\},\{u_{\Gamma}\}\right)$ may yield a value outside the allowed bounds for $\Delta_{\ell}$. This problem can be addressed by imposing appropriate restrictions on $\left\{\mathcal{G}_{12;\ell}\right\}$. Finally, $\left\{\mathcal{G}_{12;\ell}\right\}$ makes the physical interpretation of the total energy expression somewhat obscure.

In this paper, the aforementioned shortcomings are resolved using the qubit parametrization (see Section III). The qubit parametrization has two important differences. First, the qubit parametrization employs an effective many-body density matrix $\bm{\rho},$ having dimension $2^{2N_{orb}}\times 2^{2N_{orb}}$, corresponding to a pure state of a $2N_{orb}$ qubit system. The $\bm{\rho}$ is constructed such that the density of $\bm{\rho}$ is the same as the physical local density, and it can be viewed as a function of $\{u_{\Gamma}\}$ and $\left\{\mathcal{G}_{12,\ell}\right\}$. Second, the qubit parametrization uses $\Delta_{\ell}\left(\left\{\mathcal{G}_{12;\ell}\right\},\{u_{\Gamma}\}\right)$ to solve $\left\{\mathcal{G}_{12,\ell}\right\}$ as a function of $\bm{\rho}$ and $\left\{\Delta_{\ell}\right\}$, reducing the total number of constraints per spin orbital from three to two.

VI.2 Review of the gauge constrained SCDA algorithm

In this section, we use original gauge constrained SCDA to evaluate $\hat{\rho}_{G3}$. It should be noted that there are several restrictions on the variational freedom of the SPD when using the gauge constrained SCDA. First, the kinetic projector must be diagonal in momentum space. Second, the interacting projector may not introduce off-diagonal terms at the single-particle level. These two restrictions guarantee that the local integer time self-energy and Green’s function are diagonal in the original basis. For Hamiltonians with density-density interactions and hopping parameters that are diagonal in the orbital index, such as the ones treated in this paper, the aforementioned variational restrictions do not limit the variational power of the SPD. There are two critical insights in the gauge constrained algorithm. First, the integer time self-energy only has non-trivial values within the time steps containing the interacting projector, and therefore $\bm{\mathcal{G}}$ only needs to be specified in the corresponding regions. Second, the gauge symmetry can be used to restrict the form of $\bm{\mathcal{G}}$.

We start by examining the gauge symmetry of $\hat{\rho}_{G3}=\hat{K}_{1}\hat{P}_{1}\hat{K}_{2}\hat{P}_{1}^{\dagger}\hat{K}_{1}^{\dagger}$, where the gauge transformation is given by $\hat{K}_{1}\rightarrow\hat{K}_{1}\hat{N}_{1}$, $\hat{P}_{1}\rightarrow\hat{N}_{1}^{-1}\hat{P}_{1}\hat{N}_{2}^{-1}$, and $\hat{K}_{2}\rightarrow\hat{N}_{2}\hat{K}_{2}\hat{N}_{2}^{\dagger}$, where $\hat{N}_{1}=\exp\left(\bm{\mu}_{1}\cdot\hat{\bm{n}}\right)$ and $\hat{N}_{2}=\exp\left(\bm{\mu}_{2}\cdot\hat{\bm{n}}\right)$. The gauge transformation can be parametrized by

[TABLE]

as explained in Section II.2. In the following, we assume $\left[\bm{\mu}_{i}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}\mu_{i,\ell}$ where $\mu_{i,\ell}$ is a real number, yielding

[TABLE]

Notice that the interacting projectors only act on the first and second time step, and therefore it is useful to split $\mathcal{G}_{\ell}$ into the following block structure:

[TABLE]

A similar block structure is adopted for $g_{\ell}$ and $S_{\ell}$.

The first step is to focus on the $A$ block, which is sufficient to determine the integer time self-energy. Similar to the case of $\hat{\rho}_{B2}$, we only need to specify $\mathcal{G}_{\ell;A}$, which is sufficient to determine $g_{\ell;A}$ and $S_{\ell;A}$ and therefore $S_{\ell}$. Moreover, similar to the derivation of Eq. (261), we have

[TABLE]

The gauge transformation can be used to further restrict the form of $\mathcal{G}_{\ell;A}$. Given that Eq. (266) involves the inverse of $\mathcal{G}_{\ell;A}$, it is more convenient to first use Eq. (259). Notice that $\left[g_{\ell}\right]_{11}=\left[g_{\ell}\right]_{22}$ given translation symmetry and the fact that the total particle number of a given spin-orbital $\ell$ commutes with $\hat{P}_{1}$. Therefore, Eq. (259) indicates that the diagonal elements of $g_{\ell}$ are invariant, while $\mu_{1;\ell}-\mu_{2;\ell}$ can be chosen such that $\left[g^{\prime}_{\ell}\right]_{12}=-\left[g^{\prime}_{\ell}\right]_{21}$. Notice that the interacting projectors are the same for the first and second time step, and therefore $\left[\mathcal{G}_{\ell}^{\prime}\right]_{11}=\left[\mathcal{G}_{\ell}^{\prime}\right]_{22}$ and $\left[\mathcal{G}_{\ell}^{\prime}\right]_{12}=-\left[\mathcal{G}_{\ell}^{\prime}\right]_{21}$. Now consider a gauge transformation with $\mu_{1;\ell}=\mu_{2;\ell}=\mu_{\ell}$, which still preserves $\left[\mathcal{G}^{\prime\prime}_{\ell}\right]_{11}=\left[\mathcal{G}^{\prime\prime}_{\ell}\right]_{22}$ and $\left[\mathcal{G}^{\prime\prime}_{\ell}\right]_{12}=-\left[\mathcal{G}^{\prime\prime}_{\ell}\right]_{21}$, and we can choose $\mu_{\ell}$ such that $\left[\mathcal{G}^{\prime\prime}_{\ell}\right]_{11}=1/2$. Therefore, after fully exploring the gauge symmetry of the SPD, it is sufficient to use just one parameter $\mathcal{G}_{\ell,12}$ to parametrize $\mathcal{G}_{\ell;A}$ as

[TABLE]

and correspondingly

[TABLE]

It is useful to define a new quantity

[TABLE]

which is different from the $A$ block of $S_{0,\ell}=\mathcal{G}_{\ell}^{-1}-1$, denoted as $S_{0,\ell;A}$. A similar quantity $S_{F,\ell}^{\left(A\right)}$ can also be defined as

[TABLE]

which is different from the $A$ block of $S_{F,\ell}=g_{\ell}^{-1}-1$, denoted as $S_{F,\ell;A}$. Using the results of Section VI.2.2, the integer time Dyson equation within the $A$ block may be written as

[TABLE]

and $S_{\ell;A}$ can be determined as

[TABLE]

where Eqns. (267) and (268) can be used to obtain

[TABLE]

The $S_{\ell}$ is a $3\times 3$ matrix obtained from $S_{F,\ell}=S_{0,\ell}S_{\ell}$, yielding

[TABLE]

The second step is to examine the lattice integer time Green’s function and use the self-consistency condition to resolve the $B$, $C$, and $D$ blocks. The lattice integer time Green’s function can be parametrized by the physical momentum density distribution n_{k\ell}\equiv\langle\textrm{\@text@baccent{\hat{n}}}_{k\ell}^{(3)}\rangle_{\textrm{\@text@baccent{\hat{\rho}}}_{K}} and the integer time self-energy $S_{\ell}$ by assuming $\hat{K}_{2}$ is a single Slater determinant with $\left\langle\hat{n}_{k\ell}\right\rangle_{\hat{K}_{2}}=1$ for $k\in<$ and $\left\langle\hat{n}_{k\ell}\right\rangle_{\hat{K}_{2}}=0$ for $k\in>$, where the symbols $<$ or $>$ denote the occupied and unoccupied regions of $\hat{K}_{2}$. The local integer time Green’s function for the lattice is denoted as $\left[\bm{g}^{\prime}_{i}\right]_{\ell\ell^{\prime}}=\delta_{\ell\ell^{\prime}}g^{\prime}_{\ell}$, where

[TABLE]

where $g^{\prime}_{\ell,A}$ , $g^{\prime}_{\ell,B}$, and $g^{\prime}_{\ell,C}$ are $2\times 2$, $2\times 1$, and $1\times 2$ matrices, respectively, defined as

[TABLE]

where $\mathcal{A}_{<\ell}$, $\mathcal{A}_{>\ell}$, and $\Delta_{\ell}$ are defined in Eqns. (68), (69), and (67).

Using Eqns. (275), (331), (332), and (333), the $B$, $C$, and $D$ blocks of Eq. (265) can be fully determined, completing the algorithm. Finally, the local interaction energy can then be written using the tensor product representation as

[TABLE]

while the kinetic energy is $\sum_{\ell}\int dk\epsilon_{k\ell}n_{k\ell}$. The $n_{k\ell}$ are subject to the two linear constraints given in Eqns. (255) and (256), which can be implemented in the $u$-representation using

[TABLE]

while the $n_{k\ell;0}$ are subject to $\int_{<}dk=n_{\ell}$.

It is useful to appreciate what can be gleaned from the form of the preceding equations. The $g^{\prime}_{\ell,A}$ block is reminiscent of the $\hat{\rho}_{G2}$ case where $\Delta_{\ell}$ captures the quasi-particle renormalization, while $g^{\prime}_{\ell,B}$ and $g^{\prime}_{\ell,C}$ are reminiscent of the $\hat{\rho}_{B2}$ case where $\mathcal{A}_{<\ell}$ and $\mathcal{A}_{>\ell}$ capture the super-exchange effects. This observation helps illustrate how the $\hat{\rho}_{G3}$ simultaneously captures the physics of both $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$.

In the following, we further elaborate on two points which were not fully elucidated in Ref. Cheng2023035127 . First, we provide explicit expressions for the tensor product representation in the case of $\mathcal{N}=3$. Second, we explore various relations derived using the block structure of the integer time Dyson equation.

VI.2.1 Evaluating observables under \textrm{\@text@baccent{\hat{\rho}}}_{loc} using the tensor

product representation

Consider the tensor product representation of ̱\hat{O} evaluated under \textrm{\@text@baccent{\hat{\rho}}}_{loc}. Given that the interacting projector only acts on the first and second time step, the 3-dimensional tensor representation may be reduced to a 2-dimensional tensor representation as

[TABLE]

Furthermore, the tensor contraction can be simplified as (\textrm{\@text@baccent{\hat{O}}})_{u}\cdot\left(u_{1}\otimes u_{2}\right)=u^{T}(\textrm{\@text@baccent{\hat{O}}})_{u}u with $u^{T}=u_{1}=u_{2}$. Therefore, even though the compound space of $\mathcal{N}=3$ is larger than $\mathcal{N}=2$ for a given original Hilbert space, the computational cost in the tensor representation for $\hat{\rho}_{G3}$ is similar to that of $\hat{\rho}_{G2}$ due to this dimension reduction. Similar to the case of the $\hat{\rho}_{G2}$, \left(\textrm{\@text@baccent{\hat{O}}}\right)_{u} has a direct product structure given by Eq. (38). We first specify the components of the$A$ block quantities at a given spin orbital $\ell$

[TABLE]

where $S$ indicates the symmetric part of the matrix, which is defined as

[TABLE]

which is useful given that u^{T}\left(\textrm{\@text@baccent{\hat{O}}}\right)_{u}u=u^{T}\left(\textrm{\@text@baccent{\hat{O}}}\right)_{u,S}u. Similarly, the $w$-representation can be defined using $u=Vw$, where $V=\otimes_{\ell}V_{\ell}$ with

[TABLE]

such that \left(\textrm{\@text@baccent{\hat{1}}}\right)_{w;\ell} is the identity matrix and \left(\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\dagger\left(1\right)}\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\left(1\right)}\right)_{w;\ell,S} is diagonal, resulting in

[TABLE]

and correspondingly, we have

[TABLE]

Finally, we provide explicit expressions for \left(\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\dagger\left(3\right)}\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\left(3\right)}\right)_{w;\ell}, with the components given by

[TABLE]

which can be used to compute the local interaction energy.

VI.2.2 Block structure of the integer time Dyson equation

Here we derive various useful equations using the block form of the integer time Dyson equation. To make our discussion general, we assume that the integer time Green’s functions are not diagonal in the orbital index, resulting in the following block matrix equation

[TABLE]

where the blocks for $S_{F}$ and $S_{0}$ can be explicitly expressed in terms of the blocks of $g$ and $\mathcal{G}$ using the inverse formula for a $2\times 2$ block matrix, resulting in the following relations

[TABLE]

Using Eqns. (318), (319), (320), and (321), we obtain

[TABLE]

Similarly, using Eqns. (318), (319), and (322), we obtain

[TABLE]

Using Eqns. (326), (327), (324), and (325), we verify that $\mathcal{G}_{A}$ and $g_{A}$ can be used to determine $S_{A}$ as

[TABLE]

Using Eqns. (327) and (328), we have

[TABLE]

Finally, using Eq. (328), we have

[TABLE]

which can be used to solve for $\mathcal{G}_{A}$. Additionally, it is useful to explicitly write expressions for $\mathcal{G}_{B}$ , $\mathcal{G}_{C}$, and $\mathcal{G}_{D}$ as

[TABLE]

VI.3 Derivation of the qubit energy form

In this section, we derive the qubit energy form, given in Eq. (62), in several steps. First, in Section VI.3.1, we introduce a polar representation for the $A$ block. Second, in Section VI.3.2, we introduce $\bm{\rho}$ as a unitary transformation of $ww^{T}$, and we solve the self-consistency condition for the $A$ block using $n_{\ell}$, $\xi_{\ell}$, and $\Delta_{\ell}$. Third, in Section VI.3.3, we resolve the self-consistency of the $B$, $C$, and $D$ blocks using the extra information provided by $\mathcal{A}_{<\ell}$ and $\mathcal{A}_{>\ell}$.

VI.3.1 Polar Representation of the A block

We first review some mathematical properties of the matrix group with following form

[TABLE]

where $c>0$ and $\phi\in\left[0,2\pi\right]$ and the group multiplication is

[TABLE]

There is an isomorphism for this matrix group to $\mathbb{R}^{+}\times S^{1}$ with $\mathcal{S}\left(c,\phi\right)\rightarrow\left(c,\phi\right)$, with the group product taken as

[TABLE]

Interpreting $\mathcal{S}\left(c,\phi\right)$ as the integer time self-energy, the corresponding integer time Green’s function is given as

[TABLE]

and it can be seen that

[TABLE]

It is also useful to express $c$ and $\phi$ in terms of $g_{11}$ and $g_{12}$ as

[TABLE]

We now use the preceding results to study the $A$ block of the integer time Green’s function. Given that $g_{\ell,A}$ and $\mathcal{G}_{\ell,A}$ have the form given by Eq. (337), the $S_{F,\ell}^{\left(A\right)}$ and $S_{0,\ell}^{\left(A\right)}$ also have the form of Eq. (334), and therefore $S_{\ell,A}$ also has the form of Eq. (334). We can then use the polar representation to express the following quantities

[TABLE]

allowing for the integer time Dyson equation to be recast as

[TABLE]

which yields Eq. (77). Using $g_{11}\rightarrow n_{\ell}$, $g_{12}\rightarrow g_{\ell,12}$, $c\rightarrow c_{\ell},$ and $\phi\rightarrow\phi_{\ell}$ in Eqns. (340) and (341), we obtain Eqns. (75) and (76). Using $g_{11}\rightarrow 1/2$ , $g_{12}\rightarrow\mathcal{G}_{\ell,12}$, $c\rightarrow 1$, and $\phi\rightarrow\phi_{\ell,0}$, we obtain

[TABLE]

which yields Eq. (78) by assuming $\phi_{\ell,0}\in\left[-\pi/2,0\right]$. Using Eq. (344), we can rewrite Eq. (275) in terms of $c_{\ell}$, $\theta_{\ell}$ and $n_{k\ell}$ as

[TABLE]

The matrix elements in $w$ representation are given by

[TABLE]

VI.3.2 Resolving the self-consistency in the $A$-block

A key ingredient of the qubit parametrization is the introduction of the qubit representation, defined by

[TABLE]

where $\mathcal{U}$ is a unitary matrix defined as $\mathcal{U}=\mathcal{U}_{1}\otimes\dots\otimes\mathcal{U}_{2N_{orb}}$ and

[TABLE]

where $\psi_{\ell}$ is a parameter that will be determined by demanding that $\bm{\rho}$ yields the physical density. This should be contrasted to the $w$ representation, where the local density of $ww^{T}$ is generally different from the physical density.

We begin by deriving Eq. (73), which determines $\theta_{\ell}$ for a given $n_{\ell},$ $\xi_{\ell}$, and $\Delta_{\ell}$. Using Eq. (362), Eqns. (356)-(361) are transformed to

[TABLE]

where

[TABLE]

Using the self-consistency condition for the 1,1 and 1,2 entries, we have

[TABLE]

Eqns. (373) and (374) reduce to two linear equations in $p_{\ell}$ and $q_{\ell}$ as

[TABLE]

and $p_{\ell}$ and $q_{\ell}$ can be determined from Eqns. (88) and (89). Using Eqns. (371) and (372), we obtain

[TABLE]

Both $\cot\left(\phi_{\ell}\right)$ and $\cot\left(\phi_{\ell,0}\right)$ can be expressed in terms of $\cot\left(\theta_{\ell}\right)$ as

[TABLE]

Substituting Eqns. (88) and (89) into Eq. (380) and using Eq. (383), we obtain a sixth order equation in $\cot\left(\theta_{\ell}\right)$ which can be factored into the following form

[TABLE]

Notice that the first factor in Eq. (384) is positive, implying that the second factor is zero, and $\theta_{\ell}$ may be obtained as

[TABLE]

which yields Eq. (73). Notice that $\left|\cos\left(\theta_{\ell}\right)\right|\leq 1$, yielding a constraint on $\xi_{\ell}$ as

[TABLE]

VI.3.3 Resolving the $B$, $C$, and $D$ blocks

We begin by deriving a subtle symmetry between $\mathcal{G}_{B}$ and $\mathcal{G}_{C}$, described by Eq. (83) and Eq. (84), which can be rewritten as

[TABLE]

where the tilde is defined by the following rules. For a $2\times 1$ matrix $m$, we have $\left[\tilde{m}\right]_{21}=\left[m\right]_{11}$ and $\left[\tilde{m}\right]_{11}=\left[m\right]_{21}$. For a $1\times 2$ matrix $m$, we have $\left[\tilde{m}\right]_{12}=\left[m\right]_{11}$ and $\left[\tilde{m}\right]_{11}=\left[m\right]_{12}$. For a $2\times 2$ matrix $m$, we have $\left[\tilde{m}\right]_{11}=\left[m\right]_{22}$, $\left[\tilde{m}\right]_{22}=\left[m\right]_{11}$, $\left[\tilde{m}\right]_{12}=\left[m\right]_{21}$, and $\left[\tilde{m}\right]_{21}=\left[m\right]_{12}$. In order to prove Eq. (387), we use Eq. (327) to obtain

[TABLE]

Given that $\mathcal{G}_{B}=\mathcal{G}_{A}g_{A}^{-1}g_{B}$ (see Eq. (331)) and that $S_{A}$ commutes with $g_{A}$ (see Section VI.3.2), Eq. (387) can be proven by verifying that

[TABLE]

Using Eq. (331), we obtain

[TABLE]

To simplify these two equations, we introduce $\mathcal{I}_{\ell}$, $\mathcal{J}_{\ell}$, $\mathcal{A}^{\prime}_{<,\ell}$, and $\mathcal{A}^{\prime}_{>,\ell}$ (defined in Eqns. (79)-(82)), which yields Eqns. (83) and (84). To compute $\mathcal{G}_{\ell,33}$, Eqns. (333), (326), and (329) can be used to obtain

[TABLE]

Introducing $i_{\ell}$ and $j_{\ell}$, defined in Eqns. (85) and (86), we can simplify Eq. (397) to (87). The matrix elements of \left(\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\dagger\left(3\right)}\textrm{\@text@baccent{\hat{a}}}_{\ell}^{\left(3\right)}\right)_{w,\ell} can then be determined as

[TABLE]

Eqns. (362) and (364) may be used to evaluate

[TABLE]

and we define

[TABLE]

where $f_{\ell,0}$, $f_{\ell,x}$, and $f_{\ell,z}$ are given in Eqns. (90), (92), and (91). The derivations of Eqns. (71)-(92) are now complete.

VI.4 Examining the qubit energy form in special cases

In this section, we showcase the qubit energy form for the special case of half-filled orbitals. Additionally, we examine how the qubit energy form for $\hat{\rho}_{G3}$ recovers the qubit energy forms for $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$.

VI.4.1 The case of half-filled orbitals

In this section, we examine the case of half-filling with particle-hole symmetry where $n_{\ell}=1/2$ and $\mathcal{A}_{<\ell}=\mathcal{A}_{>\ell}=\frac{1}{2}A_{\ell}$. Using the general algorithm given in Eqs. (71)-(92), we provide corresponding results. Starting with $\xi_{\ell,0}=\frac{1}{2}$ and $\delta n_{\ell}=0$, we obtain

[TABLE]

and $j_{\ell}=0$, $\mathcal{G}_{\ell,33}=\frac{1}{2}$, and

[TABLE]

and $q_{\ell}=0$, $f_{\ell,0}=\frac{1}{2}$, $f_{\ell,x}=0$, and

[TABLE]

It is also useful to define

[TABLE]

yielding

[TABLE]

Eq. (428) can be plugged into Eq. (62) to obtain the total energy. This result will be applied to the multiorbital Hubbard model in Section VII.

VI.4.2 Recovering the qubit energy form for $\hat{\rho}_{G2}$

Given that $\hat{\rho}_{G2}$ is a special case of $\hat{\rho}_{G3}$, it is clear that the former can be obtained by constraining the latter. We previously demonstrated that restricting the momentum density distribution to be flat in each region and taking $\mathcal{G}_{\ell,12}=1/2$ within $\hat{\rho}_{G3}$ will recover $\hat{\rho}_{G2}$ in the case where $\{n_{kl;0}\}$ corresponds to a pure state Cheng2023035127 . Here we illustrate this fact using the qubit parametrization. We begin by enforcing $\mathcal{G}_{\ell,12}=1/2$, and Eq. (78) yields $\phi_{\ell,0}=-\pi/2$. Using Eq. (383) and $\phi_{\ell,0}=-\pi/2$, we obtain

[TABLE]

Using Eqns. (73) and (429), we obtain

[TABLE]

yielding $\mathcal{I}_{\ell}=\frac{1}{2}\left(n_{\ell}+\xi_{\ell}\right)$ and $\mathcal{J}_{\ell}=\frac{1}{2}\left(\xi_{\ell}-n_{\ell}\right)$. Using the assumption of a flat momentum density distribution, given by $n_{k\ell}|_{k\in<}=\frac{\xi_{\ell}^{2}}{n_{\ell}}+n_{\ell}$ and $n_{k\ell}|_{k\in>}=\frac{\xi_{\ell}^{2}}{n_{\ell}-1}+n_{\ell}$, which yields a quasi-particle weight of

[TABLE]

we obtain

[TABLE]

and $i_{\ell}=1$, $j_{\ell}=0$, $\mathcal{G}_{\ell,33}=\frac{1}{2}$, $p_{\ell}=-1$, $q_{\ell}=0$, $f_{\ell,0}=\frac{1}{2}$, $f_{\ell,z}=-\frac{1}{2}$, $f_{\ell,x}=0$, and

[TABLE]

yielding $\hat{n}_{eff,\ell}=\hat{n}_{\ell}$, recovering the qubit energy form obtained from $\hat{\rho}_{G2}$.

VI.4.3 Recovering the qubit energy form for $\hat{\rho}_{B2}$

Given that $\hat{\rho}_{B2}$ is a special case of $\hat{\rho}_{G3}$, it is clear that the former can be obtained by constraining the latter. However, it should be noted that in the qubit parametrization, $\hat{K}_{2}$ in $\hat{\rho}_{G3}$ is assumed to correspond to a Slater determinant, but $\hat{K}_{2}$ must be the identity to recover $\hat{\rho}_{B2}$. Nonetheless, we demonstrate that the qubit energy form of $\hat{\rho}_{G3}$ can still be constrained to recover the qubit energy form of $\hat{\rho}_{B2}$. The solution is to restrict $\bm{\rho}$ to be a diagonal matrix, implying that $\xi_{\ell}=0$, yielding $\theta_{\ell}=\frac{\pi}{2}$, $g_{\ell,12}=0$, $c_{\ell}=\frac{1}{n_{\ell}}-1$, $\phi_{\ell}=0$, $\phi_{\ell,0}=-\frac{\pi}{2}$, $\mathcal{G}_{\ell,12}=\frac{1}{2}$, $\mathcal{I}_{\ell}=\frac{n_{\ell}}{2}$, $\mathcal{J}_{\ell}=-\frac{n_{\ell}}{2}$, $\mathcal{A}^{\prime}_{<,\ell}=\mathcal{A}_{<,\ell}+\mathcal{A}_{>,\ell}=A_{\ell}$, $\mathcal{A}^{\prime}_{>,\ell}=0$, $\mathcal{G}_{\ell,13}=\mathcal{G}_{\ell,23}=\frac{A_{\ell}}{2\xi_{\ell,0}}$, $i_{\ell}=\frac{A_{\ell}^{2}}{\xi_{\ell,0}^{2}}$, $j_{\ell}=0$, $\mathcal{G}_{\ell,33}=n_{\ell}-\frac{A_{\ell}^{2}\delta n_{\ell}}{\xi_{\ell,0}^{2}}$, $p_{\ell}=-1$, $q_{\ell}=0$, f_{\ell,0}=n_{\ell}-\frac{A_{\ell}^{2}\text{\deltan}_{\ell}}{\xi_{\ell,0}^{2}}, $f_{\ell,z}=-\frac{A_{\ell}^{2}}{2\xi_{\ell,0}^{2}}$, and $f_{\ell,x}=0$. Finally, we have

[TABLE]

recovering the qubit energy form of $\hat{\rho}_{B2}$. The preceding result demonstrates that when all orbitals have $\xi_{\ell}=0$, implying that the system is a Mott insulator, the $\hat{\rho}_{G3}$ solution can be recovered by $\hat{\rho}_{B2}$, demonstrating the power of $\hat{\rho}_{B2}$ in the Mott phase.

VII Applications: multiorbital Hubbard model at half filling with particle-hole

symmetry in $d=\infty$

Here we showcase how to use the qubit energy form for $\hat{\rho}_{G3}$ to study the multiorbital Hubbard model in $d=\infty$ with density-density interactions. While the qubit energy form can be applied to arbitrary densities, here we study the case of half-filling where the local interaction energy takes a simple form. In Section VII.1, we explicitly evaluate the qubit energy form for half-filling and demonstrate how to analytically minimize over the momentum density distribution, yielding a final energy form which can straightforwardly be numerically minimized. In Section VII.2, we demonstrate that $\Delta_{\ell}$, defined in Eq. (67), is a key variable for understanding the Mott transition. For a special density-of-states, the minimization can be analytically performed, yielding an analytical relation between the ground state energy and the Hubbard $U$ via $\Delta_{\ell}$. In Section VII.3, we study the effect of the Hund coupling $J$ on the nature of the Mott transition in the multiorbital Hubbard model.

VII.1 Numerical minimization of the qubit trial energy

VII.1.1 Evaluating the qubit trial energy

In order to elucidate the results of the multiorbital Hubbard model, we first begin by considering the single orbital case with density of states $D\left(\epsilon\right)$ having particle-hole symmetry. Due to spin symmetry, we omit the spin-orbital index $\ell$ in \xi_{\ell},$$\Delta_{\ell}, and $A_{\ell}=\mathcal{A}_{<\ell}+\mathcal{A}_{>\ell}$. The qubit energy form for $\hat{\rho}_{G3}$ is given as (see Section VI.4.1 for derivation)

[TABLE]

where particle-hole symmetry implies $n\left(\epsilon\right)+n\left(-\epsilon\right)=1$, and in order to ensure that $\mathcal{F}(\Delta,A,\xi)$ is real, $\xi$ must satisfy the following constraint

[TABLE]

The effective many-body density matrix $\bm{\rho}$ can be encoded as $\bm{\rho}=|\Psi\rangle\langle\Psi|$ where $|\Psi\rangle$ can be parametrized using a single parameter $\mathcal{D}\in[0,1/2]$ as

[TABLE]

This parametrization ensures that $n_{\ell}=1/2$ and $\xi$ can be determined as

[TABLE]

The double occupancy is given as

[TABLE]

and to have a real $d(\Delta,A,\mathcal{D})$, the $h(\Delta,\mathcal{D})$ must be real which requires $\mathcal{D}$ to satisfy the following constraint

[TABLE]

Finally, the trial energy for the single orbital case is a function only of $n(\epsilon)$ and $\mathcal{D}$, given as

[TABLE]

We now proceed to the multiorbital case, with a local interaction given as

[TABLE]

where $\delta\hat{n}_{\ell}=\hat{n}_{\ell}-\frac{1}{2}$. The qubit trial energy is given as

[TABLE]

where $\bm{\rho}$ is a pure state that is restricted to $\langle\hat{n}_{\ell}\rangle_{\bm{\rho}}=\frac{1}{2}$, particle-hole symmetry implies $n_{\ell}\left(\epsilon\right)+n_{\ell}\left(-\epsilon\right)=1$, and $\xi_{\ell}$ must satisfy the following constraint

[TABLE]

The expression for $\mathcal{F}(\Delta_{\ell},A_{\ell},\xi_{\ell})$ is given in Eq. (447). The most straightforward approach to satisfying Eq. (462) is to first choose $\bm{\rho}$, and then Eq. (462) becomes a linear constraint on $n_{\ell}(\epsilon)$ (see Section III for further discussion).

VII.1.2 Numerical minimization of the qubit trial energy

We now proceed to minimizing the qubit trial energy. As before, we first focus on the single orbital model for clarity, and then consider the multiorbital case. The general numerical minimization has been described in Section III, which consists of two steps: First, the momentum density distribution is partially optimized under the constraint of $n_{\ell}$, $\Delta_{\ell}$, $\mathcal{A}_{<\ell}$, and $\mathcal{A}_{>\ell}$, or through four Lagrange multiplier $a_{<\ell}$, $a_{>\ell},$ $b_{<\ell}$, and $b_{>\ell}$. Second, one needs to minimize over the remaining $2^{2N_{orb}}+3\times 2N_{orb}-1$ independent variational parameters subjected to the inequality constraint given by Eq. (70). For the half-filled, particle-hole symmetric case, there are several simplifications. First, the optimized momentum density distribution now only depends on $\Delta_{\ell}$ and $A_{\ell}=\mathcal{A}_{<\ell}+\mathcal{A}_{>\ell}$, or just two Lagrange multipliers $a=a_{<\sigma}=-a_{>\sigma}$ and $b=b_{>\sigma}=b_{<\sigma}$, and the partially optimized momentum density distribution is given as

[TABLE]

where the Lagrange multipliers $a$ and $b$ are determined by $\Delta$ and $A$ via inverting the following relation

[TABLE]

which yields $a\left(\Delta,A\right)$ and $b\left(\Delta,A\right)$. In order to better appreciate the flexibility of $n\left(\epsilon,a,b\right)$, it is useful to examine the behavior for various choices of $a$ and $b$ (see Figure 1). For $a>0$ and $b>0$, the distribution corresponds to a Fermi liquid, given that the discontinuity of $n\left(\epsilon,a,b\right)$ at $\epsilon=0$ yields a quasi-particle weight $Z$ of

[TABLE]

which will be between zero and one. For the special case where $a\rightarrow\infty$ and $b\rightarrow\infty$ with $a/b$ remaining finite, the distribution recovers a flat distribution, which is obtained for an optimized $\hat{\rho}_{G2}$. Finally, when $a=0$ and $b>0$, the system is in the Mott phase where $Z=0$.

The second simplification for half-filling and particle-hole symmetry is that the number of independent parameters is $2^{2N_{orb}}+2N_{orb}-1$, given that $n_{\ell}=1/2$ and $\mathcal{A}_{<\ell}=\mathcal{A}_{>\ell}$. If we have further symmetry between different spin orbitals, the number of independent parameters can further be reduced. For example, for the single-orbital case with spin symmetry, we have $A_{\ell}=A$ and $\Delta_{\ell}=\Delta$, and therefore the number of independent parameters is $2^{2}+2-1-2=3$. As discussed in Section III, there are two possible strategies. The first strategy starts from $a,b$, which determines $n\left(\epsilon\right)$, and $\bm{\rho}$ is specified via $\mathcal{D}$, which will have a range given by Eq. (457). Mathematically, this qubit trial energy is given by

[TABLE]

This strategy allows the energy to be explicitly written in terms of the variational parameters, and thus it is straightforward to compute the derivatives via automatic differentiation. For the multiorbital case, the constraint between $\bm{\rho}$ and $\{\Delta_{\ell}\}$ is given by Eq. (462), which is not straightforward to implement for a given $\{\Delta_{\ell}\}$. Therefore, it is useful to pursue a second strategy.

We begin by illustrating the second strategy in the single-orbital case. This strategy begins by specifying $\bm{\rho}$ via $\mathcal{D}$, then the range of $\Delta$ can be determined from the constraint given by Eq. (451), and then $\Delta$ and $b$ can be specified, which determines $a$ via inverting Eq. (465), denoted as $a\left(\Delta,b\right)$. Mathematically, the resulting qubit trial energy is given by

[TABLE]

For the multiorbital case, the qubit trial energy is given as

[TABLE]

where $a_{\ell}(\Delta_{\ell},b_{\ell})$ and $A_{\ell}(a_{\ell},b_{\ell})$ are the multiorbital versions of $a(\Delta,b)$ and $A(a,b)$. The second strategy allows one to automatically implement all of the constraints in the multiorbital case, and the only downside is that $a_{\ell}(\Delta_{\ell},b_{\ell})$ must be numerically evaluated, though this is a trivial task.

VII.2 Understanding the Mott transition in the single orbital model

In the preceding section, the resulting qubit trial energy for the single orbital model at half-filling is either parametrized by $a$, $b$, and $\mathcal{D}$ (i.e. Eq. (467)), or $\Delta$, $b$ and $\mathcal{D}$ (i.e. Eq. (468)), which both allow for straightforward numerical minimization for a given $U$ and $D(\epsilon)$. In this section, we explore a different perspective based on the one body reduced density matrix viewpoint (see Ref. companion ), which provides a clearer understanding of the Mott transition. The first step begins with the parametrization of the qubit trial energy using $\Delta$, $A$, and $\mathcal{D}$, where the local interaction energy is $Ud(\Delta,A,\mathcal{D})$, where $d(\Delta,A,\mathcal{D})$ is defined in Eq. (455). It should be noted that $d(\Delta,A,\mathcal{D})$ has the form $1/4+A^{4}f(\Delta,\mathcal{D})$, and therefore we can optimize $d(\Delta,A,\mathcal{D})$ with respect to $\mathcal{D}$ for fixed $\Delta$ and $A$, yielding an optimized $\mathcal{D}$ as a function of $\Delta$, denoted as $\mathcal{D}\left(\Delta\right)$. Therefore, the interaction energy is purely a function of $\Delta$ and $A$, which can be viewed as a one body reduced density matrix functional, yielding a simple picture of the Mott transition. Two key points should be appreciated. First, $\mathcal{D}\left(\Delta\right)$ has a critical value of $\Delta$, denoted $\Delta_{c}$, and $\mathcal{D}\left(\Delta\right)=0$ for $\Delta\geq\Delta_{c}$, which corresponds to the Mott insulator. The second point is that $\Delta_{c}$ is independent of $D(\epsilon)$. While the total energy can be obtained by numerically minimizing over $\Delta$ and $A$, further insight can be obtained by expressing $U$ and the total energy as functions of $\Delta$.

VII.2.1 Qubit trial energy in terms of $\Delta$ and $A$

As outlined in Section VII.1.2, the qubit trial energy can be parametrized as

[TABLE]

For simplicity, we study the case where $U\geq 0$. We now consider how to minimize the total energy over $\mathcal{D}$, which amounts to minimizing $d(\Delta,A,\mathcal{D})$ over $\mathcal{D}$. In order to visualize the minimization, we plot $d$ as a function of $\mathcal{D}$ for a given $\Delta$ and $A$. Given that $A$ will not influence the minimization over $\mathcal{D}$, we choose $A=\sqrt{\left(2-4\Delta\right)\Delta}$ from a flat momentum density distribution (see Figure 2). For a given $\Delta$ curve, there are four distinct sets of points denoted as $\hat{\rho}_{G2}$, $\hat{\rho}_{B2}$, $\hat{\rho}_{G3}$, and $max$, where the $\hat{\rho}_{G3}$ points provide the optimized values for $\mathcal{D}$, the $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$ points provide the $\mathcal{D}$ values for $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$, respectively, and a $max$ point provides the maximum value for $\mathcal{D}$ given by $\mathcal{D}_{-}$ in Eq. (458). For small values of $\Delta$, the optimized value of $\mathcal{D}$ is nonzero, and $\mathcal{D}$ monotonically decreases with increasing $\Delta$. For $\Delta$ larger than some critical value, the optimized value for $\mathcal{D}$ is zero. Having obtained a graphical understanding of this function, we proceed to mathematically minimize the $d$ over $\mathcal{D}$ for a given $\Delta$ and $A$, which is a constrained minimization given that $\mathcal{D}\in[0,\mathcal{D}_{-}]$. Solving $\partial d/\partial\mathcal{D}=0$ yields

[TABLE]

When $\mathcal{D}^{\star}(\Delta)\in[0,\mathcal{D}_{-}]$, then $\mathcal{D}^{\star}(\Delta)$ yields the optimized value for $\mathcal{D}$, and otherwise the optimized value is given by the minimum value for the boundary points. It is useful to solve $\mathcal{D}^{\star}(\Delta)=0$, yielding

[TABLE]

Therefore, the optimized value for $\mathcal{D}$ is

[TABLE]

having two distinct regimes as a function of $\Delta$ (see inset of Figure 2). Finally, the physical double occupancy can be written as a function of $\Delta$ and $A$ as

[TABLE]

Therefore, the total trial energy can be written purely in terms of $\Delta$ and $A$, given as

[TABLE]

In order to find the ground state for a given $D(\epsilon)$ and $U$, it is necessary to minimize over $\Delta$ and $A$, which cannot be performed analytically in general. However, it is important to appreciate that the optimized value of $\Delta$ is sufficient to determine if the system is in the Mott phase. For $\Delta=0$, the ansatz corresponds to the Hartree-Fock wave function, while the maximum value of $\Delta=1/4$ corresponds to a collection of isolated atoms, and the Mott transition occurs when $\Delta=\Delta_{c}$, before the system becomes a collection of atoms. We now verify that the system is indeed a Mott insulator for $\Delta>\Delta_{c}$. First, the local interaction energy is independent of $\Delta$ for $\Delta>\Delta_{c}$ (see Eq. (476)), dictating that $a=0$ for the optimized $\Delta$. Second, Eq. (466) dictates that the quasiparticle weight is zero when $a=0$, implying a Mott insulating state. Alternatively, when $\Delta<\Delta_{c}$, the ansatz describes a metallic phase.

Following the results of Sections VI.4.2 and VI.4.3, we illustrate how the qubit energy form for $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$ can be recovered from $\hat{\rho}_{G3}$ when properly restricting the variational parameters. We begin with $\hat{\rho}_{B2}$, which can be viewed as a continuation of the Mott phase for $\hat{\rho}_{G3}$, where $\mathcal{D}=0$ and $d$ is given by $\frac{1}{4}-4A^{4}$ (see Figure 2). Alternatively, the $\hat{\rho}_{G2}$ is characterized by a flat momentum density distribution determined by $\Delta$ where $d$ and $\mathcal{D}$ are identical, given as

[TABLE]

It is well known that for $\hat{\rho}_{G2}$ the Mott transition occurs when the system becomes a collection of isolated atoms, meaning that the transition happens for $\Delta=1/4$.

VII.2.2 Solution for the two-peak density of states

In Section VII.2.1, we demonstrated that the interaction energy can be written analytically in terms of $\Delta$ and $A$, but the kinetic energy must be numerically determined in terms of $\Delta$ and $A$. While the latter is a trivial numerical problem, it still precludes a completely analytic solution for the energy in terms of $U$. Therefore, we introduce the two-peak density of states, where the kinetic energy can be analytically evaluated in terms of $\Delta$ and $A$, allowing for an analytical relation between the total energy and the Hubbard $U/|K_{0}|$ in all three ansatz, where $K_{0}<0$ is the non-interacting kinetic energy for the lattice. Specifically, the two-peak density of states is given by

[TABLE]

While this density of states is somewhat unphysical, the results retain all of the qualitative features of the Bethe lattice, in addition to being quantitatively very similar assuming the same $K_{0}$ (see Figure 3). As expected, the $\hat{\rho}_{G2}$ solution is identical for both density of states, given that $\hat{\rho}_{G2}$ only depends on $K_{0}$, independent of the details of $D(\epsilon)$. Interestingly, $\hat{\rho}_{G3}$ has a very similar critical value of $U/|K_{0}|$ for both density of states, with small differences differences in the double occupancy for a given $U/|K_{0}|$.

We now proceed to analytically evaluate the total trial energy as a function of $\Delta$ for the two-peak model. Using Eqns. (480), (449), and (448), the kinetic energy and $A$ are determined as $K=\left(1-4\Delta\right)K_{0}$ and $A=\sqrt{\left(2-4\Delta\right)\Delta}$, yielding the total trial energy in terms of the single variational parameter $\Delta$, given as

[TABLE]

where

[TABLE]

and $d_{1}(\Delta)$ is defined in Eq. (97). The $d(\Delta)$ is plotted in Figure 4$a$. While for a given $U$ the energy cannot be analytically minimized over $\Delta$, it is straightforward to find $U$ for a given value of $\Delta$ that satisfies $dE(\Delta)/d\Delta=0$, given as

[TABLE]

and $U/|K_{0}|$ is plotted as a function of $\Delta$ in Figure 4$b$. The critical value of $U$ and double occupancy at the Mott transition, where $\Delta=\Delta_{c}$, are given as

[TABLE]

These values are the same from both the metallic and insulating sides of the Mott transition, which can be seen in Figure 4, confirming that the transition is continuous.

Corresponding equations for $d$ and $U/|K_{0}|$ as functions of $\Delta$ can be obtained for $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$ by substituting the corresponding $d(\Delta)$ relations into Eq. (483) (see Figure 4 for plots). For $\hat{\rho}_{G2}$, we have

[TABLE]

which recovers the Gutzwiller approximation, and the Mott transition occurs at $\Delta=1/4$, where $U_{c,G2}/\left|K_{0}\right|=8$ and $d_{c,G2}=0$. For $\hat{\rho}_{B2}$, there are some subtleties to consider. For a given $U$, there are several candidate values of $\Delta$: the value given by the saddle point Eq. (483), which yields

[TABLE]

and the boundary values of $\Delta=0$ and $\Delta=1/2$. The total energy must be used to evaluate these candidate values of $\Delta$ and select the global minimum. It should be noted that $\Delta=0$ recovers the Hartree-Fock solution. There are two critical points to consider: the local stability of $\hat{\rho}_{B2}$ and a transition from a saddle point solution to the Hartree-Fock solution. The local stability of the $\hat{\rho}_{B2}$ is determined by the minimal value of $U$ in Eq. 489, which is given by

[TABLE]

For any $U>U_{c,B2}$, there exists a locally stable $\hat{\rho}_{B2}$ solution. However, one should compare the energy of the saddle-point solution to the Hartree-Fock solution, which yields another critical value of $U_{c,B2}^{\prime}/|K_{0}|=3.375$. For $U>U_{c,B2}^{\prime}$, the saddle point solution is the global minimum, while for $U<U_{c,B2}^{\prime}$ the Hartree-Fock solution is the global minimum. We summarize all the critical points for $\hat{\rho}_{G3}$, $\hat{\rho}_{G2}$, and $\hat{\rho}_{B2}$ in Table 1.

VII.2.3 Solution for a general density of states

In Section VII.2.2, we demonstrated that the qubit trial energy can be written solely in terms of $\Delta$ for the case of a two-peak density of states. Here we extend this strategy for a general density of states, allowing for the optimized value of $A$ to be evaluated as a function of $\Delta$, which is completely independent of $U$. Therefore, this approach is particularly useful for analyzing the Mott transition. We begin by rewriting the qubit trial energy from Eq. (478) using the result of Eq. (482), yielding

[TABLE]

where $\tilde{O}\left(\Delta\right)$ is given in Eq. (96). We proceed by constructing the saddle point equations of $\Delta$ and $A$ for a given $U$, yielding the following two equations

[TABLE]

where $a$ and $b$ are the Lagrange multipliers from Eq. (463). A practical approach for solving the two preceding equations for a given $U$ is to express $A$ and $\Delta$ in terms of $a$ and $b$, denoted as $A(a,b)$ and $\Delta(a,b)$, and then solve for $a$ and $b$. However, we take an alternative approach, as our goal is to determine the optimized value of $A$ for a given $\Delta$. Therefore, we proceed by dividing Eq. (494) by Eq. (495), which yields

[TABLE]

Moreover, $a$ and $b$ are required to yield the given value of $\Delta$, such that

[TABLE]

Simultaneously solving Eqns. (496) and (497) yields $a$ and $b$ as functions of $\Delta$, and therefore all quantities that depend on $n(\epsilon)$, including $A$ and $K$, are now functions purely of $\Delta$. Subsequently, the double occupancy and $U$ can be determined as a function of $\Delta$ as

[TABLE]

Alternatively, the $U$ can also be expressed as

[TABLE]

which will recover Eq. (483) when evaluating the two peak model. For the case of the Bethe lattice, we plot $A$, $a$, and $b$ as functions of $\Delta$ in Figure 5, and the critical values of all quantities are listed in Table 2.

VII.3 Understanding the effect of Hund’s coupling in the multiorbital Hubbard

model

Here we generalize the treatment from Section VII.2 to the multiorbital case including the Hund coupling $J$. The key difference is that in the multiorbital case, one cannot analytically minimize over the local variational parameters, though these parameters can easily be numerically minimized as a function of $\Delta$ for a given $J/U$. The remaining procedure closely follows the single orbital case.

Consider the multiorbital Hubbard model defined in Eq. (94). The qubit trial energy can be written as

[TABLE]

where $\mathcal{F}\left(\Delta,A,\xi\right)$ is defined in Eq. (447), $\xi(\bm{\rho})=\frac{1}{2}\langle\hat{\sigma}_{\ell}^{x}\rangle_{\bm{\rho}}$, and

[TABLE]

where $\hat{O}_{i}$ are defined in Eq. (94). Having defined the qubit trial energy, we proceed to obtain the solutions for the two-peak density of states and the Bethe lattice using the $\Delta$ parametrization. We first rewrite the local interaction energy as

[TABLE]

where

[TABLE]

which is obtained from $\mathcal{F}\left(\Delta,A,\xi\right)/A^{2}$. Given the form of Eq. (505), the optimized $\bm{\rho}$ will only depend on $\Delta$, motivating the definition of the following function

[TABLE]

A convenient way to generate this function is to perform a two stage minimization. First, we perform a constrained minimization with the restriction $\xi\left(\bm{\rho}\right)=\xi$ on $\bm{\rho}$. Second, we minimize the expression over $\xi$. In the first stage, given that $\tilde{\mathcal{F}}^{2}\left(\Delta,\xi\left(\bm{\rho}\right)\right)$ is fixed, we only need to minimize $\langle\hat{O}\rangle_{\bm{\rho}}$, which can be mathematically expressed as

[TABLE]

To efficiently generate $\mathcal{O}\left(\xi\right)$ , we can introduce a Lagrange multiplier $\lambda$ and determine the ground state for $\hat{\mathcal{H}}=\hat{O}-\lambda\sum_{\alpha\sigma}\hat{\xi}_{\alpha\sigma}$, yielding the optimal $\bm{\rho}$ and corresponding $\xi$ and $\langle\hat{O}\rangle_{\bm{\rho}}$ for a given $\lambda$. One can then perform such calculations over a grid of $\lambda$, and then spline the relationship between $\langle\hat{O}\rangle_{\bm{\rho}}$ and $\xi$. A plot of $\mathcal{O}(\xi)$ is provided in Figure 2 of Ref. companion . Finally, the partially optimized local energy can be written as

[TABLE]

In order to visualize $\tilde{O}\left(\Delta\right)$ and $\xi$ as a function of $\Delta$, we plot these quantities for $J/U=0,0.05,0.25$ and $N_{orb}=2,3,5,7$ (see Figure 6). For $J/U=0$, one can clearly observe that $\xi$ continuously goes to zero, while there is a discontinuity for $J/U=0.05$. For $J/U=0.25$, $\xi$ discontinuously goes to zero for $N_{orb}=2,3$ and continuously goes to zero for $N_{orb}=5,7$. It should be noted that $\tilde{O}\left(\Delta\right)$ can be applied to solve an arbitrary particle-hole symmetric $D(\epsilon)$, and therefore captures the essence of the Mott transition for a given $J/U$.

VII.3.1 Understanding the non-analytic behavior of $\tilde{O}\left(\Delta\right)$

via a Taylor series

In Section VII.2, we demonstrated that there is non-analyticity in $\tilde{O}\left(\Delta\right)$ at $\Delta=\Delta_{c}$. Here we provide a Taylor series analysis to explain how the non-analyticity emerges. Equation (510) indicates that $\tilde{O}\left(\Delta\right)$ is the minimum of $\tilde{\mathcal{F}}^{2}\left(\Delta,\xi\right)\mathcal{O}(\xi)$ within the range $\xi\in\left[0,\frac{1}{2}-\Delta\right]$, and it is convenient to study the quantity

[TABLE]

Finding the minimum of $\mathcal{L}\left(\Delta,\xi\right)$ will yield the minimum of $\tilde{\mathcal{F}}^{2}\left(\Delta,\xi\right)\mathcal{O}(\xi)$ given that $\tilde{\mathcal{F}}\left(\Delta,0\right)=4$ and $\mathcal{O}\left(0\right)=-\frac{1}{4}N_{orb}\left(1+\left(N_{orb}-1\right)J/U\right)$.

We begin by Taylor series expanding $\mathcal{L}\left(\Delta,\xi\right)$ to sixth order in $\xi$ about $\xi=0$, and the second order coefficient in $\xi$ is expanded in $\Delta$ about $\Delta_{c,I}$ such that it is zero for $\Delta=\Delta_{c,I}$, yielding

[TABLE]

where $c_{2}$, $c_{4}$, $c_{6}$, and $\Delta_{c,I}$ are constants for a given $N_{orb}$ and $J/U$. Given that the optimized $\xi$ goes to zero with increasing $\Delta$ (see Figure 6), this requires $c_{2}>0$. Furthermore, we take $c_{6}>0$, though there will be cases where $c_{6}$ is negative and a higher order expansion is necessary. Therefore, we only need to understand how the sign of $c_{4}$ influences the non-analyticity in $\tilde{O}\left(\Delta\right)$. For $c_{4}>0$, when $\Delta>\Delta_{c;I}$, the minimum of $\mathcal{L}$ is given by $\mathcal{L}=-1$ with $\xi=0$, while for $\Delta<\Delta_{c;I}$, the minimum of $\mathcal{L}$ is obtained with

[TABLE]

Notice that $\xi^{2}$ continuously increases from [math] when $\Delta$ decreases from $\Delta_{c;I}$. Therefore the critical value of $\Delta$ is given by $\Delta_{c}=\Delta_{c;I}$. Moreover, using

[TABLE]

we find that $\frac{d\tilde{O}\left(\Delta\right)}{d\Delta}$ is continuous and only has a kink at $\Delta=\Delta_{c}$. For $c_{4}<0$, when $\Delta>\Delta_{c;I}$ there is a local minimum at $\xi=0$, while for $\Delta<\Delta_{c,I}+\frac{c_{4}^{2}}{3c_{2}c_{6}}$ there is a local minimum given by Eq. (513). The two saddle points need to be compared to obtain the global minimum, which yields $\Delta_{c}=\Delta_{c,I}+\frac{c_{4}^{2}}{4c_{2}c_{6}}.$ Therefore, $\xi$ jumps from zero to a finite value when $\Delta$ decreases from $\Delta_{c}$, implying that $\frac{d\tilde{O}\left(\Delta\right)}{d\Delta}$ is discontinuous at $\Delta=\Delta_{c}$. It should be noted that for $c_{4}>0$, $\Delta_{c}=\Delta_{c;I}$ is exact, while for $c_{4}<0$, the expression $\Delta_{c}=\Delta_{c,I}+\frac{c_{4}^{2}}{4c_{2}c_{6}}$ is an approximation, and in this case one should use the exact form of $\mathcal{L}\left(\Delta,\xi\right)$ to determine $\Delta_{c}$ if precision is needed.

We now proceed to analytically compute the expansion for two cases: $J/U=0$ with $N_{orb}\geq 1$ , and $J/U>0$ with $N_{orb}=2$ . We begin by expanding

[TABLE]

where

[TABLE]

It should be noted that $a_{2}$, $a_{4}$, and $a_{6}$ all monotonically decrease with increasing $\Delta$ for $\Delta\in\left[0,1/4\right]$. It is straightforward to show that $a_{2}=0$ when $\Delta=1/4$ and $a_{4}=0$ when $\Delta=\frac{1}{4}\left(2-\sqrt{\frac{1}{3}\left(2+\sqrt{7}\right)}\right)\approx 0.188895$ and $a_{6}=0$ when $\Delta=0.155281$. The remaining task is to compute $-\frac{\mathcal{O}\left(\xi\right)}{\mathcal{O}\left(0\right)}$, which we separately consider for the two aforementioned cases.

For $J/U=0$ and $N_{orb}\geq 1$, symmetry can be used to parametrize $\bm{\rho}=|\Psi\rangle\langle\Psi|$ with $|\Psi\rangle=\sum_{\Gamma}\sqrt{p_{\left|N_{\Gamma}-N_{orb}\right|}}|\Gamma\rangle$, where $p_{i}$ is a variational parameter with $i=0,\dots,N_{orb}$ and $N_{\Gamma}=\langle\Gamma|\sum_{\ell=0}^{2N_{orb}}\hat{n}_{\ell}|\Gamma\rangle$ counts the number of electrons in state $\Gamma$. It is convenient to reparametrize $p_{i}$ as

[TABLE]

such that $\text{Tr}\left(\bm{\rho}\right)=\sum_{i=0}^{N_{orb}}x_{i}^{2}$, where $(m)_{n}=\frac{\Gamma(m+n)}{\Gamma(m)}=m(m+1)...(m+n-1)$ is the Pochhammer symbol. When taking an expectation value of an operator $\hat{A}$ in the qubit space, it is convenient to use a matrix representation $\left[\hat{A}\right]_{ij}$, where $i$ and $j$ take values from $0,\dots,N_{orb}$, such that

[TABLE]

The non-zero entries of the effective matrices for $\hat{O}$ and $\hat{\sigma}_{\ell}^{x}$ are given as

[TABLE]

[TABLE]

Perturbation theory can then be used to obtain

[TABLE]

where

[TABLE]

The expansion coefficients of $\mathcal{L}$ are obtained as

[TABLE]

Given that $c_{4}>0$, we have $\Delta_{c}=\Delta_{c;I}$ and $\frac{d\tilde{O}\left(\Delta\right)}{d\Delta}$ has a kink at $\Delta=\Delta_{c}$.

We now discuss the case of $J/U>0$ with $N_{orb}=2$. Using perturbation theory, we find

[TABLE]

where

[TABLE]

where $r=J/U$. Therefore, we can compute the coefficients in the expansion of $\mathcal{L}$ as

[TABLE]

We can see that $c_{4}<0$, and therefore $\Delta_{c}>\Delta_{c,I}$, indicating that there is a discontinuity in $\frac{d\tilde{O}\left(\Delta\right)}{d\Delta}$ for $\Delta=\Delta_{c}$.

Finally, we numerically explore the cases of $N_{orb}>2$, where $\Delta_{c;I}$ and $c_{2}$ are identical to the case of $N_{orb}=2$. To be concrete, we consider $J/U=0.25$, and we numerically compute the expansion coefficients by fitting $\mathcal{O}(\xi)$ to a sixth order polynomial. For $N_{orb}=3$, we have $c_{4}\approx-2.9$ and $c_{6}\approx-147$. Notice that in this case $c_{6}<0$, and thus an expansion beyond sixth order is necessary, though plotting $\mathcal{L}\left(\Delta,\xi\right)$ indicates that $\Delta_{c}>\Delta_{c,I}$, yielding a discontinuity in $\frac{d\tilde{O}\left(\Delta\right)}{d\Delta}$ for $\Delta=\Delta_{c}$ (not shown). For $N_{orb}=4$, $c_{4}\approx-0.63$ and $c_{6}\approx 39$. For $N_{orb}=5$, $c_{4}\approx 0.44$ and $c_{6}\approx 35$. For $N_{orb}=6$, $c_{4}\approx 1.1$ and $c_{6}\approx 33$. For $N_{orb}=7$, $c_{4}\approx 1.5$ and $c_{6}\approx 34$. It should be noted that for $N_{orb}\geq 5$, we have $\Delta_{c}=\Delta_{c;I}$ and $\frac{d\tilde{O}\left(\Delta\right)}{d\Delta}$ has a kink at $\Delta=\Delta_{c}$.

VII.3.2 Solution for the two-peak density of states

We now consider the two peak density of states

[TABLE]

where $K_{0}$ is the total non-interacting kinetic energy per site. For the two-peak density of states, $A=\sqrt{\Delta(2-4\Delta)}$, and the qubit trial energy can be written purely in terms of $\Delta$ as

[TABLE]

Following the single orbital case, $\Delta$ can be used to determine $U$ from $dE\left(\Delta\right)/d\Delta=0,$ which yields

[TABLE]

thus providing a succinct solution parametrized by $\Delta$. The relation $U(\Delta)$ can be used to determine the nature of the Mott transition from $\widetilde{O}\left(\Delta\right)$, given that this quantity will allow any observable to be expressed in terms of $U$. We plot $\frac{dO_{tp}}{d\Delta}$ versus $\Delta$ for various $N_{orb}$ and $J/U$, which demonstrates three types of non-analytical scenarios for $O_{tp}\left(\Delta\right)$ (see Fig. 7, panel $a$). First, for all $J/U=0$, the $\frac{dO_{tp}}{d\Delta}$ is continuous with a positive slope and has a kink at $\Delta=\Delta_{c}$. Second, for $J/U>0$ and small $N_{orb}$, the $\frac{dO_{tp}}{d\Delta}$ is discontinuous at $\Delta=\Delta_{c}$ with a negative slope for $\Delta_{c}^{-}$. Third, for $J/U=0.25$ and $N_{orb}=7$, the $\frac{dO_{tp}}{d\Delta}$ is continuous and has a kink at $\Delta=\Delta_{c}$, and the slope is negative for $\Delta_{c}^{-}$. The $\frac{dO_{tp}}{d\Delta}$ can now be used to determine $U(\Delta)$, which yields the order of the Mott transition (see Fig. 7, panel $b$). First, for all $J/U=0$, $U$ increases monotonically and continuously with $\Delta$, with a kink at $\Delta=\Delta_{c}$, and therefore there are no metastable regions and the Mott transition is continuous. Second, for $J/U>0$, there is an unstable region in the metal phase where $\frac{dU}{d\Delta}<0$, and the total energy can be used to determine the transition between the metal and insulating phase, corresponding to a horizontal line, and the transition is first-order.

In summary, we have demonstrated that the nature of the Mott transition for the two peak density of states is determined purely by $\widetilde{O}\left(\Delta\right)$, and below we demonstrate that the Bethe lattice has the same behavior. Therefore, it appears that $\widetilde{O}\left(\Delta\right)$ is the essence of what determines the nature of the Mott transition in $d=\infty$.

VII.3.3 Solution for a general density of states

We now execute a similar strategy for a general density of states. Using Eq. (509), the qubit trial energy can be written as

[TABLE]

The saddle point equations are given as

[TABLE]

Given that $\Delta$ and $A$ are functions of $a$ and $b$, one can solve $a$ and $b$ from Eqns. (544) and (545) for a given $U$, and then determine all physical quantities. We now demonstrate that $\frac{d\widetilde{O}\left(\Delta\right)}{d\Delta}=0$ indicates that the system is in the Mott phase. Given that the quasiparticle weight is given as $Z=a/\sqrt{a^{2}+b^{2}}$ (see Eq. 466), and that $A>0$ for finite $U$ (see Eq. 464), the only scenario where $Z=0$ is when $\frac{d\widetilde{O}\left(\Delta\right)}{d\Delta}=0$. Therefore, when $\Delta<\Delta_{c}$, the system is metallic, while $\Delta>\Delta_{c}$ the system is insulating. For a given $U$, one must minimize over $\Delta$ in order to determine nature of the ground state.

An alternate approach is to parametrize the solution in terms of $\Delta$. Equations (544) and (545), in addition to the constraint on $a$ and $b$ for a given $\Delta$, yield

[TABLE]

For a given $\Delta$, $a$ and $b$ can be determined from Eqns. (546) and (547), $A$ can be determined from $A\left(a,b\right)$, and $U$ can be determined using Eq. (545) as

[TABLE]

allowing for the evaluation of the total energy.

We now consider the Bethe lattice in $d=\infty$ for $N_{orb}=2,3,5,7$. Recall that for a given $J/U$, the $\tilde{O}\left(\Delta\right)$ yields a $\Delta_{c}$ which divides the metallic and insulating states, where $\Delta>\Delta_{c}$ indicates an insulating phase. For the case of a continuous transition, the $U_{c}$ will be determined by $\Delta_{c}$, while for a first-order transition, one must explicitly determine the $U_{c}$ where the insulating and metallic states cross in energy. The algorithm is executed by evaluating $A$, $U$, the interaction energy, and the total energy as functions of $\Delta$. We begin by plotting the $U/t$ as a function of $\Delta$ (see Figure 8, panel $a$). For $J/U=0$, $U$ is a monotonic and continuous function of $\Delta$, implying a continuous phase transition at $\Delta_{c}$, which can be identified as a kink. Alternatively, for $J/U=0.05$, the $U$ is not a monotonic nor a continuous function of $\Delta$, implying that the there are regions of phase coexistence and unstable regions. The metallic curve exists for $\Delta<\Delta_{c}$, and the solution is only stable for $\Delta<\Delta_{c;1}$, where $\Delta_{c;1}$ is determined from $\frac{dU}{d\Delta}=0$, and therefore the metallic solution is only stable for $U<U(\Delta_{c;1})$. The insulating curve exists for $\Delta>\Delta_{c}$ without any unstable regions, and therefore the insulating phase exists for $U>U(\Delta_{c}+0^{+})$. Given that $U(\Delta_{c}+0^{+})<U(\Delta_{c;1})$, there exists a region of coexistence for the metallic and insulating solutions, and the total energy dictates the lowest energy solution. We now proceed to present the total energy and the interaction energy as functions of $U/t$ for various $J/U$ and $N_{orb}$ (see Figure 9), and the results are very similar to the two-peak case. Consistent with previous Gutzwiller Bunemann19974011 , slave boson Hasegawa19971391 , and DMFT Ono2003035119 studies, the Mott transition is continuous for $J/U=0$ and first-order for $J/U>0$.

Finally, we evaluate $U_{c}$ for $J/U=0$, which is explicitly given by $U_{c}=\frac{b_{c}}{8A_{c}^{3}}$, where $b_{c}$ is determined from $\Delta(0,b_{c})=\Delta_{c}$ and $A_{c}=A(0,b_{c})$, where $\Delta(a,b)$ and $A(a,b)$ are defined in Eqns. (465) and (464), respectively. The critical values for the Bethe lattice are listed in Table 3. For the case of large $N_{orb}$, we have

[TABLE]

consistent with our numerical results in Ref. Cheng2023035127 . For a systematic exploration of how $U_{c}$ depends on $J/U$ and $N_{orb}$, see Fig. 5 of Ref. companion .

VIII Summary and Conclusions

We begin by providing a high level overview of the variational discrete action theory (VDAT), such that the developments of the present work can be properly understood. VDAT is a variational approach to the many-body body problem that consists of two main components: a variational ansatz for the many-body wave function or density-matrix, known as the sequential product density matrix (SPD), and a formalism for evaluating expectation values under the SPD, known as the discrete action theory Cheng2021195138 ; Cheng2021206402 . The SPD has a natural mechanism to trade off between efficiency and accuracy, where the integer $\mathcal{N}$ monotonically increases the variational power of the SPD and guarantees the ability to recover the ground state solution for $\mathcal{N}\rightarrow\infty$. Moreover, there are two distinct types of SPD which satisfy the properties of a many-body density matrix, denoted as G-type and B-type. The G-type $\mathcal{N}=1$, $\mathcal{N}=2$, and $\mathcal{N}=3$ SPD encapsulate the Hartree-Fock wave function, the Gutzwiller wave function, and the Gutzwiller-Baeriswyl wave function, respectively. The key breakthrough using VDAT was the demonstration that the SPD can be exactly evaluated for multiorbital Hubbard models in $d=\infty$. We demonstrated that the G-type $\mathcal{N}=3$ SPD accurately solves the Anderson impurity model on a ring Cheng2021206402 , the single band Hubbard model over all parameter space Cheng2021206402 , the two orbital Hubbard model including a crystal field and the full rotationally invariant Hund’s coupling Cheng2022205129 , and the $SU(2N_{orb})$ Hubbard model for $N_{orb}\leq 8$ Cheng2023035127 . Moreover, we demonstrated that the computational cost of solving a G-type $\mathcal{N}=3$ SPD is comparable to a G-type $\mathcal{N}=2$ SPD, meaning that VDAT can provide a sufficiently accurate solution at a cost not far beyond the Gutzwiller approximation. The success of VDAT at $\mathcal{N}=3$ motivated a search for the best possible algorithm for executing calculations using a G-type $\mathcal{N}=3$ SPD Cheng2023035127 , which is essential for detailed exploration of the multiorbital Hubbard model and merging VDAT with realistic electronic structure methods.

The VDAT algorithm in $d=\infty$ consists of two steps: the exact evaluation of the SPD via the self-consistent canonical discrete action theory (SCDA) and the optimization of the energy with respect to the variational parameters. The SCDA requires the numerical solution of a set of self-consistency conditions, and therefore can be inconvenient when minimizing over the variational parameters. For the case of a G-type $\mathcal{N}=3$ SPD with certain restrictions (see Sections I, VI.1, and VI.2 for further details), the SCDA self-consistency condition can be automatically satisfied, which we refer to as the gauge constrained SCDA algorithm Cheng2023035127 . In the present work, we introduce the so-called qubit parametrization of the gauge constrained SCDA algorithm, which is mathematically equivalent to the original gauge constrained SCDA algorithm. The qubit parametrization offers several key improvements. The qubit parametrization analytically resolves some constraints over the variational parameters, thus reducing the number of variational parameters by one per spin orbital. Additionally, the variational parameters are physically intuitive and facilitate a deeper understanding of how the SPD captures Mott and Hund physics. Therefore, the qubit parametrization achieves the long sought goal of resolving the shortcomings of the Gutzwiller approximation while maintaining the computational simplicity and physical appeal.

The variational parameters of the qubit parametrization consist of the momentum density distribution, the non-interacting reference momentum density distribution, and the pure state of a qubit system with a dimension of the local Hilbert space. The qubit system naturally arises from reparametrizing the variational parameters of the interacting projector, and the renormalized correlations within the qubit space yield the physical local correlations. The variational parameters are restricted by two constraints per spin orbital, requiring that the local density computed from the momentum density distribution is the same as that computed from the non-interacting reference momentum density distribution, and the same as the density computed from the qubit system. The qubit trial energy has a very intuitive form: the kinetic energy is determined by the momentum density distributions, while the local interaction energy is the expectation value of an effective Hamiltonian within the qubit system. Interestingly, the effective Hamiltonian has the same form as the local interacting Hamiltonian where the local density operator is substituted by an effective density operator of the qubit system. The effective density operator for a given spin-orbital $\ell$ depends on five parameters: the density $n_{\ell}$, the magnetization in the $x$ direction for the $\ell$-th qubit, denoted $\xi_{\ell}$, and three quantities determined from the momentum density distribution, denoted as $\Delta_{\ell}$, $\mathcal{A}_{<\ell}$, and $\mathcal{A}_{>\ell}$. The quantity $\Delta_{\ell}$ characterizes the number of electrons promoted across the reference Fermi surface, while $\mathcal{A}_{<\ell}$ and $\mathcal{A}_{>\ell}$ characterize the momentum density distribution below and above the reference Fermi surface, respectively. The quantity $\xi_{\ell}$ is an important variable which differentiates between a zero and non-zero quasiparticle weight when all variables are fully optimized, where $\xi_{\ell}=0$ indicates zero quasiparticle weight for spin-orbital $\ell$. The main computational cost for evaluating the ansatz is dictated by computing expectation values within the qubit system.

While evaluating the trial qubit energy ansatz is a straightforward task, optimizing over all the variational parameters remains nontrivial. In general, the qubit trial energy can be partially optimized over the momentum density distribution by introducing four Lagrange multipliers per spin-orbital, replacing the continuous momentum density distribution with four variables. For a system with $2N_{orb}$ spin-orbitals per site, there remain $2^{2N_{orb}}+3\times 2N_{orb}-1$ variational parameters which must be optimized in general. However, this number may be greatly reduced by symmetry, and it is likely possible to compress these variables into a smaller number of parameters without a serious loss of fidelity. For the special case of half-filled orbitals with particle hole symmetry, we demonstrate that one can efficiently minimize over all variational parameters, with a computational cost proportional to computing the ground state of a Hamiltonian defined within the qubit system.

In order to demonstrate the power of the qubit parametrization, we studied the ground state properties of the multiorbital Hubbard model at half-filling with particle-hole symmetry for various $J/U$ and $N_{orb}=2-7$. For a given $J/U$, the majority of the energy minimization can be encapsulated into the computation of a single variable function $\tilde{O}(\Delta)$, which can then be used to obtain the solution at a negligible cost for an arbitrary $U$ and density-of-states. The entire function $\tilde{O}(\Delta)$ is evaluated by solving a collection of qubit systems, which has a relatively small computational cost. For example, for $N_{orb}=7$, the ground state for a given qubit system can be solved in several seconds on a typical single desktop computer core, and taking on the order of 100 samples, the entire function $\tilde{O}(\Delta)$ can be accurately obtained on the order of hundreds of seconds. The extreme computational efficiency of the qubit parametrization in this case allows one to easily map out all of parameter space, which is not possible with DMFT given the lack of efficient impurity solvers for the zero temperature multiorbital problem. We find that for $J/U=0$, the Mott transition is continuous, while it is first-order for $J/U>0$, consistent with previous Gutzwiller Bunemann19974011 , slave boson Hasegawa19971391 , and DMFT Ono2003035119 studies.

While the key result of this paper is formulating the qubit parametrization of the gauge constrained SCDA algorithm for a G-type $\mathcal{N}=3$ SPD, we also demonstrate that the qubit parametrization can be applied to the G-type and B-type $\mathcal{N}=2$ SPD. Moreover, we demonstrate that properly restricting the variational parameters of the qubit trial energy for the G-type $\mathcal{N}=3$ SPD can recover the corresponding qubit trial energy for the G-type and B-type $\mathcal{N}=2$ SPD. Interestingly, the qubit trial energy for the G-type $\mathcal{N}=2$ SPD has an identical form to the slave spin mean-field theory (see Appendix B), and thus the $\mathcal{N}=3$ qubit trial energy may provide insights for proceeding beyond mean-field theory in the slave spin formalism.

The qubit parametrization of the gauge constrained SCDA algorithm at $\mathcal{N}=3$ is likely the optimal form when evaluating an SPD with a kinetic projector that is diagonal in both the momentum and spin-orbital indices and an interacting projector that consists of diagonal Hubbard operators. For the Hamiltonians treated in the present study, which have density-density interactions and hopping parameters that are diagonal in the spin-orbital index, the aforementioned restrictions on the SPD do not limit the variational power. When solving a general Hamiltonian which includes the full rotationally invariant form of the Hund exchange or non-diagonal hopping terms, the qubit parametrization can still be applied and it will still yield an upper bound on the energy in $d=\infty$, but it will not contain the full variational power of $\mathcal{N}=3$. Ongoing research is addressing how to generalize the qubit parametrization to handle an arbitrary $\mathcal{N}=3$ G-type SPD, with aspirations of completely superseding our general decoupled minimization algorithm for $\mathcal{N}=3$ Cheng2022205129 .

IX Acknowledgments

This work was a supported by a RISE-LDRD grant from Columbia University and Brookhaven National Laboratory. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

Appendix A One-body reduced density matrix functional for the Hubbard model

An important application of the qubit energy form for $\hat{\rho}_{G3}$ is the construction of a one body reduced density matrix functional (1RDMF) for the multi-orbital Hubbard model, which is the focus of our companion manuscript companion . Here we derive the corresponding results for $\hat{\rho}_{G2}$ and $\hat{\rho}_{B2}$. Additionally, we evaluate existing 1RDMF’s from the literature for the single band Hubbard model at half-filling in $d=\infty$.

A.1 The 1RDMF’s for the multi-orbital Hubbard model from $\hat{\rho}_{G2}$

and $\hat{\rho}_{B2}$

We begin by presenting the 1RDMF from $\hat{\rho}_{B2}$ using the qubit parametrization. The interaction energy is given as

[TABLE]

where $n_{\ell}=\int dkn_{k\ell}$, the effective density operator is $\hat{n}_{eff,\ell}=n_{\ell}+\frac{A_{\ell}^{2}}{\xi_{\ell,0}^{2}}\left(\hat{n}_{\ell}-n_{\ell}\right)$ with $A_{\ell}=\int dk\sqrt{n_{k\ell}\left(1-n_{k\ell}\right)}$ and $\xi_{\ell,0}=\sqrt{n_{\ell}(1-n_{\ell})}$, and both $\bm{\rho}$ and $H_{loc}\left(\left\{\hat{n}_{eff,\ell}\right\}\right)$ are diagonal in the Pauli-Z basis of the qubit system. For the case of half-filling, the minimization yields

[TABLE]

where $\mathcal{O}\left(0\right)=-\frac{1}{4}N_{orb}\left(1+\left(N_{orb}-1\right)J/U\right)$.

We now present the 1RDMF for $\hat{\rho}_{G2}$. Unlike $\hat{\rho}_{G3}$ or $\hat{\rho}_{B2}$, the variational parameters do not explicitly contain $n_{k\ell},$ but instead $n_{k\ell}-n_{\ell}=\xi_{\ell}^{2}/\xi_{\ell,0}^{2}\left(n_{k\ell,0}-n_{\ell}\right)$, where $n_{k\ell,0}\in\left[0,1\right]$. Therefore, for a given $\left\{n_{k\ell}\right\}$, we have $\left|\xi_{\ell}/\xi_{\ell,0}\right|\geq\sqrt{Z_{\ell,min}}$, where $Z_{\ell,min}=\max\left(\frac{n_{\ell,max}-n_{\ell}}{1-n_{\ell}},\frac{n_{\ell}-n_{\ell,min}}{n_{\ell}}\right)$, where $n_{\ell,max}$ and $n_{\ell,min}$ are the maximum and minimum values of $\left\{n_{k\ell}\right\}$ for a given $\ell$. It should be noted that the gauge symmetry can be used to restrict $\xi_{\ell}\geq 0$. The interaction energy can then be written as

[TABLE]

where $\xi_{\ell;min}=\xi_{\ell;0}\sqrt{Z_{min}}$. Notice that the interaction energy will decrease with decreasing $\xi_{\ell}$, and therefore the restriction in Eq. (555) can be replaced with $\frac{1}{2}\langle\hat{\sigma}_{\ell}^{x}\rangle_{\bm{\rho}}=\xi_{\ell;min}$ . For the case of half-filling, we have

[TABLE]

where $\mathcal{O}(\xi)$ is defined in Eq. (508).

Using the above interaction energy functionals will yield results identical to the corresponding VDAT results, such as the ones provided in Fig. 3.

A.2 Results using published 1RDMF’s for the single-orbital Hubbard model

in $d=\infty$

We are not aware of any applications of existing 1RDMF’s to the Hubbard model in $d=\infty$, though there are various studies of Hubbard clusters Kamil2016085141 ; Mitxelena2017425602 and $d=1$ Lopezsandoval20001764 ; Mitxelena20201701 and $d=2$ Hubbard models Saubanere2016045102 ; Mitxelena2020064108 . It is important to emphasize that $d=\infty$ is the most relevant test of local electronic correlations, and is emblematic of typical three dimensional strongly correlated materials in that $d=\infty$ hosts a standard Fermi liquid in the metal phase and exhibits a Mott transition at a finite value of $U$. While $d=1$ is exactly solvable via the Bethe ansatz, the Mott transition occurs at an infinitesimal $U$ Lieb19681445 . Alternatively, $d=2$ is a testbed for the most advanced and expensive computational approaches, and the properties at half-filling and low temperatures are still actively studied Tanaka2019205133 ; Schafer2021011058 ; Chatzieleftheriou2024236504 ; Geng2025115143 . While it is still very interesting to compare total energies from 1RDMF’s in $d=1$ and $d=2$, if the goal is to determine whether or not a 1RDMF can describe the Mott transition, then it is most critical to first benchmark in $d=\infty$.

In this section, we discuss how to use various 1RDMF’s to solve the one band Hubbard model at half-filling, including the MBB Muller1984446 , CA Csanyi20007348 , CGA Csanyi2002032510 , power Sharma2008201103 , PNOF5 Piris2011164102 , PNOF7 Piris2017063002 , and dimer Saubanere2016045102 functionals. The interaction energy is given by $Ud$, where $d$ is the double occupancy and is a functional of $\left\{n_{k\sigma}\right\}$. If the interaction Hamiltonian is rewritten as $\frac{U}{2}\sum_{i\sigma\sigma^{\prime}}\hat{a}_{i\sigma}^{\dagger}\hat{a}_{i\sigma^{\prime}}^{\dagger}\hat{a}_{i\sigma^{\prime}}\hat{a}_{i\sigma}$, the interaction energy for the MBB, CA, CGA, and power functionals can be viewed as a modification of the Hartree-Fock (HF) energy where the Fock term is altered Mitxelena2017425602 . We begin by writing the local interaction as $E_{int}=\frac{U}{2}\sum_{\sigma\sigma^{\prime}}\left\langle\hat{a}_{i\sigma}^{\dagger}\hat{a}_{i\sigma^{\prime}}^{\dagger}\hat{a}_{i\sigma^{\prime}}\hat{a}_{i\sigma}\right\rangle$, which can be written in momentum space as

[TABLE]

where $L$ is the number of $k$-points. The MBB, CA, CGA, and power functionals are all defined using the following approximation

[TABLE]

where $n_{k}$ is defined through $\left\langle\hat{a}_{k\sigma}^{\dagger}\hat{a}_{k^{\prime}\sigma^{\prime}}\right\rangle=n_{k}\delta_{kk^{\prime}}\delta_{\sigma\sigma^{\prime}}$. The $\mathcal{F}\left(n_{i},n_{j}\right)$ and $E_{int}$ are given for each functional in Table 4, and our values of $\mathcal{F}\left(n_{i},n_{j}\right)$ are identical to those provided in Ref. Mitxelena2017425602 . It is worth noting that the MBB, CA, and CGA recover the atomic limit where $E_{int}=0$ when $n_{k}=1/2$, while the power functional does not if $\alpha\neq 1/2$.

Now we discuss how to minimize the total energy over $\{n_{k}\}$. Given that MBB is a special case of the power functional, we first consider the power functional. Given that the interaction energy for the power functional is fixed for a given $\int dkn_{k}^{\alpha}$ and $\int dkn_{k}$, the kinetic energy can be minimized by introducing two Lagrange multipliers $a$ and $b$, yielding a target functional $K=\int dk\left(\epsilon_{k}-a\right)n_{k}-b\int dkn_{k}^{\alpha}$ to be minimized, and the partially optimized $n_{k}$ are obtained as

[TABLE]

The ground state energy for a given $U$ is then obtained by optimizing the total energy over $a$ and $b$.

We now consider the CA and CGA, which can be viewed in a unified way by writing the interaction energy functional as

[TABLE]

where $\theta=1$ corresponds to CA and $\theta=2$ corresponds to CGA. We can similarly introduce two Lagrange multipliers, yielding a target functional

[TABLE]

where $\theta\geq 1$ and the partially optimized value of $n_{k}$ is given as

[TABLE]

For a given $U$, one can then optimize the Lagrange multipliers to obtain the ground-state energy and other observables. Plots of the double occupancy as a function of $U/t$ are provided in Ref. companion .

An important drawback of the MBB, power, CA, and CGA functionals is that the anti-symmetry of the two-body reduced density matrix (2RDM) is violated when $\mathcal{F}\left(n_{i},n_{j}\right)$ deviates from the HF value Mitxelena2017425602 . The Piris natural orbital functionals (PNOF) Piris2013620 resolve this issue by approximating the form of the 2RDM $D_{ij,kl}^{\alpha\beta}=\frac{1}{2}\langle\hat{a}_{i\alpha}^{\dagger}\hat{a}_{j\beta}^{\dagger}\hat{a}_{l\beta}\hat{a}_{k\alpha}\rangle$ as Mitxelena2017425602

[TABLE]

where $\Delta_{ij}$ and $\Pi_{ik}$ are explicitly defined in PNOF5 and PNOF7. In the previous calculations Mitxelena2017425602 ; Piris20131298 , it was found that the optimized natural orbital often breaks symmetry in order to further minimize the energy, and translational symmetry breaking in the Hubbard model allows PNOF5 and PNOF7 to recover the correct atomic limit Mitxelena2017425602 . In our context, we require translational symmetry to be respected, so we will examine PNOF5 and PNOF7 in this case. PNOF5 only includes intra-pair correlation (i.e. $\Pi_{ik}=0$ when $i,k$ belong to different pairs), and therefore the deviation of the summation in Eq. (557) from the corresponding Hartree-Fock value scales like $L$ in the thermodynamic limit while there is a prefactor of $1/L^{2}$, and thus PNOF5 will yield the same energy form as HF in the thermodynamic limit. PNOF7 accounts for interpair correlation, which will provide a thermodynamic contribution, and for half-filling we find $E_{int}=\frac{1}{4}-\frac{1}{2}A^{2}$, which is similar to the CA but with different prefactor for $A^{2}$.

Finally, we discuss the interaction energy of the dimer functional at half-filling, given as

[TABLE]

where $K=2\int dk\epsilon_{k}n_{k}$ and $K_{0}$ is the non-interacting kinetic energy. Minimizing the total energy over $\{n_{k}\}$ yields the double occupancy as a function of $U$ as

[TABLE]

Appendix B Equivalence between the Qubit energy form for $\hat{\rho}_{G2}$

and the Slave Spin Mean Field Theory

Here we prove that the slave spin mean-field theory (SSMF) is identical to the qubit energy form for $\rho_{G2}$. We begin by describing the SSMF using the conventions in our work, as there are trivial differences. We associate the spin state $|\uparrow\rangle$ with the fermionic state $|0\rangle$, while the standard SSMF associates $|\downarrow\rangle$ with $|0\rangle$. To restore the enlarged Hilbert space of the SSMF back to the physical space, we require that $\hat{n}_{i\ell}=\frac{1-\hat{\sigma}_{i\ell}^{z}}{2}$ holds for every local site $i$ and spin-orbital index $\ell$. Then, the local Hamiltonian at site $i$ can be written as $H_{loc}\left(\left\{\hat{n}_{i\ell}\right\}\right)=H_{loc}\left(\left\{\frac{1-\hat{\sigma}_{i\ell}^{z}}{2}\right\}\right)$, while the hopping term can be written as $\hat{a}_{i\ell}^{\dagger}\hat{a}_{j\ell}=\hat{f}_{i\ell}^{\dagger}\hat{O}_{i\ell}^{\dagger}\hat{O}_{j\ell}\hat{f}_{j\ell}$, where $\hat{O}_{i\ell}=\hat{S}_{i\ell}^{+}+c_{i\ell}\hat{S}_{i\ell}^{-}$ , $\hat{S}_{i\ell}^{\pm}=\frac{\hat{\sigma}_{i\ell}^{x}\pm\hat{\sigma}_{i\ell}^{y}}{2}$, $c_{i\ell}$ is an arbitrary complex number, and $\hat{f}_{j\ell}$ is a fermionic annihilation operator. While the preceding transformations are exact for any $c_{i\ell}$, an appropriate choice of $c_{i\ell}$ is indeed critical when making mean-field approximations. In Ref. De'medici2005205124 , which only addressed half-filling, the $c_{i\ell}=1$ and $\hat{O}_{i\ell}=\hat{\sigma}_{i\ell}^{x}$, while later work generalized this result Hassan2010035106 , choosing $c_{i\ell}$ such that the non-interacting limit can always be correctly recovered within the mean field approximation. We follow the latter choice. Finally, the Hamiltonian in the slave spin representation is given as

[TABLE]

with the operator constraint $\hat{f}_{i\ell}^{\dagger}\hat{f}_{i\ell}=\frac{1-\hat{\sigma}_{i\ell}^{z}}{2}$ for every site and spin-orbital. There are two ways to derive the SSMF. The first approach introduces a mean-field decoupling for the hopping term and treats the operator constraint in the mean field level, yielding decoupled Hamiltonians for electrons and spins that must be solved self-consistently. The second approach Georgescu2017165135 ; Maurya2021425603 ; Maurya2022055602 ; Crispino2023155149 uses the variational principle, assuming a trial wave-function $|\Psi\rangle=|\Psi_{f}\rangle\otimes|\Psi_{s}\rangle$ which is a direct product of a fermionic part $|\Psi_{f}\rangle$ and a spin part $|\Psi_{s}\rangle$. Furthermore, the $|\Psi_{f}\rangle$ is assumed to be Slater determinant and the $|\Psi_{s}\rangle$ is assumed to be a direct product in real space, such that $|\Psi_{s}\rangle=\otimes_{i}|\Psi_{s;i}\rangle$. Therefore, the trial energy of Eq. (567) is given as

[TABLE]

subject to the constraint

[TABLE]

Given translational symmetry and the correspondence $\bm{\rho}\leftrightarrow|\Psi_{s;j}\rangle\langle\Psi_{s;j}|$ and $n_{k\ell;0}\leftrightarrow\langle\Psi_{f}|f_{k\ell}^{\dagger}f_{k\ell}|\Psi_{f}\rangle$, we can prove that trial energy in Eq. (568) is identical to the qubit trial energy for $\hat{\rho}_{G2}$ given in Eq. (50) and the constraint in Eq. (569) is identical to Eq. (54). The key step is to evaluate

[TABLE]

where gauge symmetry allows $\langle\hat{\sigma}_{\ell}^{y}\rangle_{\bm{\rho}}=0$. To ensure $\langle\Psi_{s;j}|\hat{O}_{j\ell}|\Psi_{s;j}\rangle=1$ in the non-interacting limit, we must choose $c_{\ell}$ such that $\left(1+c_{\ell}\right)\xi_{\ell;0}=1$, where $\xi_{\ell;0}=\sqrt{n_{\ell}\left(1-n_{\ell}\right)}$ is the value for $\xi_{\ell}=\langle\frac{1}{2}\hat{\sigma}_{\ell}^{x}\rangle_{\bm{\rho}}$ in the non-interacting limit. Plugging $c_{\ell}$ into Eq. (570), we have $\langle\Psi_{s;j}|\hat{O}_{j\ell}|\Psi_{s;j}\rangle=\xi_{\ell}/\xi_{\ell;0}$, completing the proof.

The preceding proof demonstrates the equivalence of the qubit energy form for $\hat{\rho}_{G2}$ and the trial energy for the SSMF. Now we demonstrate how to obtain a Hamiltonian form for the SSMF using the saddle point equations for the trial energy. The total energy under a fixed $\{n_{\ell}\}$ can be minimized by introducing Lagrange multipliers for the electron and spin systems, yielding

[TABLE]

Taking the derivative with respect to $\bm{\rho}$ yields the mean-field Hamiltonian for the spin system as

[TABLE]

where $\hat{n}_{\ell}=\frac{1}{2}\left(1-\hat{\sigma}_{\ell}^{z}\right)$ , $\hat{\xi}_{\ell}=\frac{1}{2}\hat{\sigma}_{\ell}^{x}$, $h_{\ell}=\frac{1}{L}\frac{\xi_{\ell}}{n_{\ell}\left(1-n_{\ell}\right)}\sum_{k}\epsilon_{k\ell}n_{k\ell,0}$, $C$ is a constant that does not influence the results, and we assume $\sum_{k}\epsilon_{k\ell}=0$. Now consider the derivative of the electron components as

[TABLE]

where $Z_{\ell}=\frac{\xi_{\ell}^{2}}{n_{\ell}\left(1-n_{\ell}\right)}$, which can be connected to a non-interacting fermionic Hamiltonian as

[TABLE]

The ground states for Eqns. (572) and (574) will yield an updated $\bm{\rho}$ and $n_{k\ell}$, yielding new Hamiltonians for Eqns. (572) and (574), and this procedure is iterated until self-consistency is achieved.

Appendix C The Central Point Expansion (CPE)

The central point expansion (CPE) can be viewed as an approach to evaluate both a G-type and B-type SPD at $\mathcal{N}=2$ by expanding about a reference SPD referred to as the central point. The CPE can be applied in arbitrary dimensions, and it is formally exact if all orders are summed, though it has only ever been applied at first-order. Interestingly, for the G-type $\mathcal{N}=2$ SPD, the first-order CPE yields the same result as the Gutzwiller approximation (GA) and the SCDA in any dimension, and the derivation offers an alternative perspective from the GA and the SCDA. For the B-type $\mathcal{N}=2$ SPD, the first-order CPE yields the same result as the SCDA in $d=\infty$, while providing a different approximation in finite dimensions. The CPE was originally developed in the context of the off-shell effective energy theory (OET) Cheng2020081105 , where the CPE was renormalized using both weak coupling and strong coupling perturbation theory to ensure the correct limiting behavior, yielding excellent results for the Hubbard model in $d=1,2,\infty$.

C.1 The CPE for $\hat{\rho}_{G2}$

In this section, we use the first-order CPE to evaluate $\hat{\rho}_{G2}$, which will be shown to be equivalent to the GA and the SCDA. The first-order CPE can be motivated by the fact that the GA relation between $n_{k\ell}-n_{\ell}$ and $n_{k\ell;0}-n_{\ell}$ is a linear form (see Eq. (111)) given by

[TABLE]

which is valid for an arbitrary $n_{k\ell;0}\in\left[0,1\right]$ with a constraint $\left(1/L\right)\sum_{k}n_{k\ell;0}=n_{\ell}$, where $L$ is the number of sites in the lattice. This remarkably simple relation motivates the use of $n_{k\ell;0}-n_{\ell}$ as the expansion parameters, where $\mathcal{Z}_{\ell}$ can be determined using a $\hat{\rho}_{G2}$ with $n_{k\ell;0}$ slightly deviating from the uniform distribution $n_{\ell}$ at first order. Alternatively, the relation between $n_{k\ell;0}$ and $\gamma_{k\ell}$ is highly nonlinear, suggesting that a first-order approximation in $\gamma_{k\ell}$ cannot be applied to the Gutzwiller wavefunction.

The CPE begins by choosing an expansion point for $\hat{\rho}_{G2}=\hat{P}_{1}\hat{K}_{2}\hat{P}_{1}$ as $\hat{\rho}_{G2}^{\star}=\hat{P}_{1}\hat{K}_{2}^{\star}\hat{P}_{1}$, where $\hat{K}_{2}^{\star}=\exp\left(\sum_{k\ell}\gamma_{\ell}^{\star}\hat{n}_{k\ell}\right)$ and $\gamma_{\ell}^{\star}$ is chosen to reproduce the non-interacting local density $n_{\ell;0}\equiv\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2}}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2}^{\star}}$ through $n_{\ell;0}=1/\left(1+\exp\left(-\gamma_{\ell}^{\star}\right)\right)$. Considering a kinetic projector that deviates slightly from $\hat{K}_{2}^{\star}$ as $\hat{K}_{2}=\exp\left(\sum_{k\ell}\gamma_{k\ell}\hat{n}_{k\ell}\right)$, where $\gamma_{k\ell}=\gamma_{\ell}^{\star}+\delta\gamma_{k\ell}$ , we compute the response of the expectation value $\langle\hat{O}\rangle_{\hat{\rho}_{G2}}$ to $\gamma_{k\ell}$ to the first order about the central point, given as

[TABLE]

where we have utilized that fact that $\hat{K}_{2}^{\star}$ commutes with $\hat{P}_{1}$ Cheng2020081105 , $\sqrt{\hat{\rho}_{G2}^{\star}}=\hat{P}_{1}\sqrt{\hat{K}_{2}^{\star}}=\sqrt{\hat{K}_{2}^{\star}}\hat{P_{1}}$, and the notation $\langle\hat{A};\hat{B}\rangle_{\hat{\rho}}$ is defined as

[TABLE]

The response coefficient of $\langle\hat{O}\rangle_{\hat{\rho}_{G2}}$ to $\gamma_{k\ell}$ at a general $\hat{\rho}_{G2}$ can be conveniently expressed via the integer time correlation function in the compound space as

[TABLE]

While Eq. (580) is difficult to evaluate in general, the CPE circumvents this problem by observing that Eq. (580) becomes trivial to evaluate at the central point given that $\hat{\rho}_{G2}^{\star}$ and $\sqrt{\hat{\rho}_{G2}^{\star}}$ are direct product states in real space. Though it is not yet clear, it will prove critical to exploit the gauge symmetry of the SPD to restrict $\hat{P}_{1}$ such that $\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{\rho}_{G2}^{\star}}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2}^{\star}}=n_{\ell;0}$. The response for the momentum density distribution $n_{k\ell}$ to $\gamma_{k^{\prime}\ell^{\prime}}$ can be obtained as

[TABLE]

where the correlation function $\langle\hat{n}_{k\ell};\hat{n}_{k^{\prime}\ell^{\prime}}\rangle_{\hat{\rho}_{G2}^{\star}}$ can be computed by transforming to real space as

[TABLE]

Utilizing the fact that $\hat{\rho}_{G2}^{\star}=\otimes_{i}\hat{\rho}_{G2;i}^{\star}$, where $\hat{\rho}_{G2;i}^{\star}$ is the local reduced density matrix for $\hat{\rho}_{G2}^{\star}$, then $\langle\hat{a}_{j_{1}\ell}^{\dagger}\hat{a}_{j_{2}\ell};\hat{a}_{j_{3}\ell^{\prime}}^{\dagger}\hat{a}_{j_{4}\ell^{\prime}}\rangle_{\hat{\rho}_{G2}^{\star}}$ is only non-zero in the following cases: (1) when $j_{1}=j_{2}$, $j_{3}=j_{4},$ and $j_{1}\neq j_{4}$, the value is $n_{\ell;0}n_{\ell^{\prime};0}$; (2) when $j_{1}=j_{4}$, $j_{2}=j_{3}$, $j_{1}\neq j_{2}$, and $\ell=\ell^{\prime}$, the value is $A_{\ell}^{2}$, where $A_{\ell}\equiv\langle\hat{a}_{j\ell}^{\dagger};\hat{a}_{j\ell^{\prime}}\rangle_{\hat{\rho}_{G2}^{\star}}=\langle\hat{a}_{j\ell};\hat{a}_{j\ell^{\prime}}^{\dagger}\rangle_{\hat{\rho}_{G2}^{\star}}$; (3) when $j_{1}=j_{2}=j_{3}=j_{4},$ and $\ell=\ell^{\prime}$, the value is $n_{\ell;0}=A_{\ell;0}^{2}-n_{\ell;0}^{2}$, where $A_{\ell;0}\equiv\langle\hat{a}_{j\ell}^{\dagger};\hat{a}_{j\ell^{\prime}}\rangle_{\hat{K}_{2}^{\star}}=\sqrt{\left(1-n_{\ell;0}\right)n_{\ell;0}}$; (4) when $j_{1}=j_{2}=j_{3}=j_{4}$, and $\ell\neq\ell^{\prime}$, the value is $\left\langle\hat{n}_{i\ell}\hat{n}_{i\ell^{\prime}}\right\rangle_{\hat{\rho}_{G2}^{\star}}$. Combining these four cases, we have

[TABLE]

where $C_{\ell\ell^{\prime}}$ is defined as

[TABLE]

Therefore, the response coefficient of $n_{k\ell}$ to $\gamma_{k^{\prime}\ell^{\prime}}$ is given as

[TABLE]

which consists of two contributions: (1) the coherent contribution $\delta_{\ell\ell^{\prime}}\delta_{kk^{\prime}}A_{\ell}^{2}$, which reflects the hopping renormalization captured in the GA, and (2) the incoherent contribution $C_{\ell\ell^{\prime}}/L$, which has no contribution to the momentum density distribution given that the constraint $\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2}}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{K}_{2}^{\star}}$ requires

[TABLE]

When applying Eq. (585) with $\hat{P}_{1}\rightarrow 1$ such that $\hat{\rho}_{G2}^{\star}\rightarrow\hat{K}_{2}^{\star}$, the response for the bare momentum density distribution $n_{k\ell;0}$ to $\gamma_{k^{\prime}\ell^{\prime}}$ is given by

[TABLE]

The next stage is to find the relation between $\delta n_{k\ell}=\left\langle\hat{n}_{k\ell}\right\rangle_{\hat{\rho}_{G2}}-\left\langle\hat{n}_{k\ell}\right\rangle_{\hat{\rho}_{G2}^{\star}}$ and $\delta n_{k\ell;0}=\left\langle\hat{n}_{k\ell}\right\rangle_{\hat{K}_{2}}-\left\langle\hat{n}_{k\ell}\right\rangle_{\hat{K}_{2}^{\star}}$. We first compute the response of $\delta n_{k\ell}$ and $\delta n_{k\ell;0}$ for a given $\delta\gamma_{k\ell}$ with the constraint given by Eq. (586), and then solve $\delta\gamma_{k\ell}$ from $\delta n_{k\ell;0}$ and express $\delta n_{k\ell}$ in terms of $\delta n_{k\ell;0}$, given as

[TABLE]

The constraint given by Eq. (586) for $\delta\gamma_{k\ell}$ naturally yields the following relation

[TABLE]

which is the constraint imposed in the GA.

Finally, to connect Eq. (588) with Eq. (107), we need to show $\mathcal{R}_{\ell}=A_{\ell}/A_{\ell;0}$, which can be accomplished by proving that $\hat{\rho}_{G2}$ and $\hat{\rho}_{G2}^{\star}$ have the same local reduced density matrix. To prove this, we need to consider the expectation of a given Hubbard operator under $\hat{\rho}_{G2}$, which can be accomplished by considering the linear response of a general diagonal Hubbard operator $\hat{X}_{i\Gamma}$ at site $i$ to $\gamma_{k\ell}$ as

[TABLE]

Given the constraint on $\gamma_{k\ell}$ (see Eq. 586), the first order change in $\langle\hat{X}_{i\Gamma}\rangle_{\hat{\rho}_{G2}}$ is zero and therefore we have

[TABLE]

Similarly, with $\hat{P}_{1}\rightarrow 1$ and $\hat{\rho}_{G2}\rightarrow\hat{K_{2}},$we have

[TABLE]

Therefore, we have recovered Eq. (102). Moreover, we have

[TABLE]

where we have used the fact that the local reduced density matrix of $\hat{\rho}_{G2}^{\star}$ is same as the local reduced density matrix of $\hat{\rho}_{G2}$, as indicated in Eq. (595), and represented as $\rho_{loc}$ in Eq. (105). Similarly, the local reduced density matrix of $\hat{K}_{2}^{\star}$ is the same as the local reduced density matrix of $\hat{K}_{2}$, as indicated in Eq. (596), and represented as $\rho_{loc;0}$ in Eq. (106).

In summary, we have demonstrated that the first order CPE is equivalent to the GA and the SCDA. The key steps in the proof include: (1) the local reduced density matrix for $\hat{\rho}_{G2}$ is the same as for $\hat{\rho}_{G2}^{\star}$, with no dependency on the details of $n_{k\ell;0}$ except its average value and (2) the momentum density distribution $n_{k\ell}$ is uniformly shrunk towards $n_{\ell}$ through $\mathcal{Z}_{\ell}$, which is uniquely determined by $\hat{\rho}_{G2}^{\star}$. Finally, it would be interesting to explore the behavior of the CPE beyond first order in finite dimensions, which provides insight beyond the GA for evaluating $\hat{\rho}_{G2}$.

C.2 The CPE for $\hat{\rho}_{B2}$

In this section, we use the first-order CPE to evaluate $\hat{\rho}_{B2}=\hat{K}_{1}\hat{P}_{1}\hat{K}_{1}$, where $\hat{P}_{1}=\exp\left(\sum_{i\Gamma}\upsilon_{\Gamma}\hat{X}_{i\Gamma}\right)$ and $\hat{K}_{1}=\exp\left(\sum_{k\ell}\gamma_{k\ell}\hat{n}_{k\ell}\right)$, demonstrating the equivalence to the SCDA in $d=\infty$. The CPE for $\hat{\rho}_{B2}$ can be viewed as a dual version of the CPE for the $\hat{\rho}_{G2}$, which has various correspondences Cheng2020081105 . First, in the CPE for the $\hat{\rho}_{G2}$, the central projector $\hat{K}_{2}$ determines the local density, and is invariant after applying the interacting projector $\hat{P}_{1}$. Correspondingly, in the CPE for $\hat{\rho}_{B2}$, the central projector $\hat{P}_{1}$ determines the local density, and is invariant after applying $\hat{K}_{1}$. Second, in the CPE for the $\hat{\rho}_{G2}$, the physical momentum density distribution is linearly related to the bare momentum density distribution, and the local reduced density matrix is independent of the details of the momentum density distribution. Correspondingly, in the CPE for the $\hat{\rho}_{B2}$, the local reduced density matrix is linearly related to the reference local reduced density matrix, and the momentum density distribution is independent of the details of the local reduced density matrix.

We proceed in evaluating $\hat{\rho}_{B2}$ via the CPE by choosing the central point for $\hat{\rho}_{B2}$ as $\hat{\rho}_{B2}^{\star}=\hat{K}_{1}\hat{P}_{1}^{\star}\hat{K}_{1}$, where $\hat{P}_{1}^{\star}=\exp\left(\sum_{i\ell}\gamma_{\ell}^{\star}\hat{n}_{i\ell}\right)$ and $\gamma_{\ell}^{\star}$ is chosen such that $n_{\ell}^{\star}\equiv\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{P}_{1}}=\left\langle\hat{n}_{i\ell}\right\rangle_{\hat{P}_{1}^{\star}}$. We can also rewrite $\hat{P}_{1}^{\star}=\exp\left(\sum_{i\Gamma}\upsilon_{\Gamma}^{\star}\hat{X}_{i\Gamma}\right)$ by expressing $\hat{n}_{i\ell}$ as a linear combination of Hubbard operators, where $\upsilon_{\Gamma}^{\star}=\sum_{\ell}\gamma_{\ell}^{\star}\Gamma(\ell)$. Considering the linear response of $\upsilon_{\Gamma}$ about $\upsilon_{\Gamma}^{\star}$, similar to Eqns. (576)-(578), we have

[TABLE]

where we similarly utilize the fact that $\hat{P}_{1}^{\star}$ commutes with $\hat{K}_{1}$ and $\sqrt{\hat{\rho}_{B2}^{\star}}=\hat{K}_{1}\sqrt{\hat{P}_{1}^{\star}}=\sqrt{\hat{P}_{1}^{\star}}\hat{K}_{1}$. We proceed by computing the response

[TABLE]

Directly computing Eq. (601) is cumbersome, and this can be avoided by decomposing $\hat{X}_{i\Gamma}$ into an alternate form using Eq. (3), resulting in

[TABLE]

where

[TABLE]

We introduce the density fluctuation operator $\delta\hat{D}_{iI}=\prod_{\ell\in I}\delta\hat{n}_{i\ell}$, where $\delta\hat{n}_{i\ell}=\hat{n}_{i\ell}-n_{\ell}^{\star}$, and $\delta\hat{D}_{I}=\hat{1}$ when $I=\left\{\right\}$, and we use $\delta\hat{D}_{iI}$ as a new basis for the CPE expansion. The diagonal Hubbard operator $\hat{X}_{i\Gamma}$ can be written as

[TABLE]

where $I$ enumerates over all subsets of $\left\{1,2,\dots,2N_{orb}\right\}$ and $\bar{I}=\left\{1,2,\dots,N_{orb}\right\}-I$ . The interacting projector can be written as $\hat{P}_{1}=\exp\left(\sum_{iI}\eta_{iI}\delta\hat{D}_{iI}\right)$, and at the central point $\eta_{i\left\{\ell\right\}}^{\star}=\gamma_{\ell}^{\star}$, with $\eta_{i\left\{\right\}}^{\star}$ absorbing the constant contribution, while $\eta_{iI}^{\star}=0$ if $\left|I\right|>1$. Similar to Eq. (600), for $\left|I\right|>0$ we have

[TABLE]

given that $\langle\delta\hat{D}_{iI}\rangle_{\hat{\rho}_{B2}^{\star}}=0$ in this case. It should be noted that $\eta_{i\{\}}$ has no effect on the expectation values so we implicitly only consider cases with $\left|I\right|>0$. The response of $\langle\delta\hat{D}_{i^{\prime}I^{\prime}}\rangle_{\hat{\rho}_{B2}}$ to $\eta_{iI}$ is given as

[TABLE]

Since $\hat{\rho}_{B2}^{\star}$ is diagonal in $\ell$, the correlation function on the right side of Eq. (607) is only non-zero when $I=I^{\prime}$, and it can be written as the product of contributions from each relevant spin-orbital as

[TABLE]

requiring the evaluation of

[TABLE]

and $\langle\hat{n}_{i\ell};\hat{n}_{i^{\prime}\ell}\rangle_{\hat{\rho}_{B2}^{\star}}$ can be evaluated similarly to Eq. (582) as

[TABLE]

Given that $\hat{\rho}_{B2}^{\star}=\otimes_{k\ell}\hat{\rho}_{B2;k\ell}^{\star}$ is a direct product state in momentum space, $\langle\hat{a}_{k_{1}\ell}^{\dagger}\hat{a}_{k_{2}\ell};\hat{a}_{k_{3}\ell}^{\dagger}\hat{a}_{k_{4}\ell}\rangle_{\hat{\rho}_{B2}^{\star}}$ is only non-zero in the following cases. (1) $k_{1}=k_{2}$, $k_{3}=k_{4}$, and $k_{1}\neq k_{4}$ results in $n_{k_{1}\ell}^{\star}n_{k_{3}\ell}^{\star}$, where $n_{k\ell}^{\star}=\langle\hat{n}_{k\ell}\rangle_{\hat{\rho}_{B2}^{\star}}$, (2) $k_{1}=k_{4}$, $k_{2}=k_{3}$, and $k_{1}\neq k_{2}$ results in $A_{k_{1}\ell}A_{k_{2}\ell}$, where $A_{k\ell}\equiv\langle\hat{a}_{k\ell}^{\dagger};\hat{a}_{k\ell}\rangle_{\hat{\rho}_{B2}^{\star}}=\langle\hat{a}_{k\ell};\hat{a}_{k\ell}^{\dagger}\rangle_{\hat{\rho}_{B2}^{\star}}$ and $A_{k\ell}=\sqrt{n_{k\ell}^{\star}\left(1-n_{k\ell}^{\star}\right)}$, (3) $k_{1}=k_{2}=k_{3}=k_{4}$ results in $n_{k\ell}^{\star}=\left(n_{k\ell}^{\star}\right)^{2}+A_{k\ell}^{2}$. In summary, we have

[TABLE]

which can be used to obtain

[TABLE]

where $L$ is the number of sites in the lattice. A constraint is imposed on $\hat{K}_{1}$ such that

[TABLE]

resulting in $\left(1/L\right)\sum_{k}n_{k\ell}^{\star}=n_{\ell}^{\star}$ and

[TABLE]

In infinite dimensions, only local contributions need to be accounted for, which is used in the following steps. Considering the response $\delta D_{iI}=\langle\delta\hat{D}_{iI}\rangle_{\hat{\rho}_{B2}}$ and $\delta D_{iI;0}=\langle\delta\hat{D}_{iI}\rangle_{\hat{P}_{1}}$ to first order in $\eta_{iI}$ and solving for $\eta_{iI}$ as a function of $\delta D_{iI;0}$, the response in infinite dimensions is

[TABLE]

where the renormalization factor is given by

[TABLE]

where $A_{k\ell}=\sqrt{n_{k\ell}^{\star}\left(1-n_{k\ell}^{\star}\right)}$ and $A_{k\ell;0}=\sqrt{n_{\ell}^{\star}\left(1-n_{\ell}^{\star}\right)}$.

Finally, we need to consider how the momentum density distribution is influenced by $\eta_{iI}$ using the response

[TABLE]

which is only non-zero when $I=\left\{\ell\right\}$. To first order, $\eta_{i\{\ell\}}=\eta_{i\{\ell\}}^{\star}$ and the contribution from $\eta_{iI}$ where $\left|I\right|>0$ yields

[TABLE]

and

[TABLE]

In conclusion, if we write the local Hamiltonian as $\hat{H}_{loc}=\sum_{iI}E_{loc;I}\delta\hat{D}_{iI}$, the total energy per site for $\hat{\rho}_{B2}$ is given as

[TABLE]

where $\delta D_{iI;0}=\langle\delta\hat{D}_{iI}\rangle_{\hat{P}_{1}}$ with $\delta D_{iI;0}=0$ for $\left|I\right|=1$, and $n_{k\ell}\equiv\langle\hat{n}_{k\ell}\rangle_{\hat{\rho}_{B2}}\in\left[0,1\right]$ is the physical momentum density distribution and is constrained by $\left(1/L\right)\sum_{k}n_{k\ell}=n_{\ell}$. The $n_{k\ell}$ can be viewed as variational parameters determined from the $\gamma_{k\ell}$ within $\hat{K}_{1}$. Combining Eq. (615) and Eq. (618), we can explicitly express $\mathcal{F}_{\ell}$ as

[TABLE]

To evaluate the Hamiltonian in the usual Hubbard operator representation for a given $X_{i\Gamma;0}=\langle\hat{X}_{i\Gamma}\rangle_{\hat{P}_{1}}$ while respecting the density constraint, we can compute $\delta D_{iI;0}$ from $X_{i\Gamma;0}$, and then use Eq. (614) to evaluate $\delta D_{iI}$ and then use Eq. (605) to compute $X_{i\Gamma}=\langle\hat{X}_{i\Gamma}\rangle_{\hat{\rho}_{B2}}.$ Finally, we remark that Eq. (614) is derived by ignoring non-local contributions, which is equivalent to taking the limit of infinite dimensions. One could straightforwardly account for non-local contributions using Eq. (613), which will be explored in the future work.