Alleviating the Sign Problem in Quantum Monte Carlo Simulations of Spin-Orbit-Coupled Multi-Orbital Hubbard Models
Aaram J. Kim Philipp Werner Roser Valentí
cond-mat.str-el



aaram.kim@kcl.ac.uk
Department of Physics, King's College London, Strand, London WC2R 2LS, UK

 
Department of Physics, University of Fribourg, Chemin du Musée 3, 1700 Fribourg, Switzerland


Institut für Theoretische Physik, Goethe-Universität Frankfurt, Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany


	We present a strategy to alleviate the sign problem
  in continuous-time quantum Monte Carlo
(CTQMC) simulations of the dynamical-mean-field-theory (DMFT) equations
for the spin-orbit-coupled multi-orbital Hubbard model.  
We first identify the combinations of rotationally invariant Hund coupling terms present in the relativistic basis
which lead to a severe sign problem.  
Exploiting the fact that the average sign in CTQMC depends on the
choice of single-particle basis, we  propose a bonding-antibonding basis V_j3/2BA which shows an improved
average sign compared to the widely used relativistic basis for most parameter
sets investigated.  We then generalize this procedure
by introducing a stochastic optimization algorithm
that exploits the space of single-particle bases
 and show that V_j3/2BA is
very close to optimal within the  parameter space investigated. 
Our findings enable more efficient DMFT simulations of materials with strong spin-orbit coupling.  


Alleviating the sign problem in quantum Monte Carlo simulations
of spin-orbit-coupled multiorbital Hubbard models
    Roser Valentí
    Received: date / Accepted: date
=================================================================================================================





§ INTRODUCTION

Spin-orbit coupling (SOC) is an essential ingredient in the study of exotic
phases in correlated electron systems,<cit.>
 such as unconventional
superconductivity in 4d transition-metal oxides,<cit.>
topological phases of matter in quantum spin-Hall insulators,<cit.>
excitonic
insulators, <cit.> and
 Kitaev-model-based insulators,<cit.>
to mention a few.
A prototypical minimal
 model that includes the interplay between spin-orbit coupling
and correlations is the relativistic  multiorbital Hubbard model.
Its nonrelativistic counterpart has been intensively investigated in the past and shows a rich phase
diagram.<cit.>
The relativistic multiorbital Hubbard model,<cit.> however, is much less understood since 
the choice of algorithms is strongly limited 
due to the extra computational complications associated with (multiple-)spin-orbit-coupled
 degrees of freedom.

One of the promising formalisms to investigate the
Hubbard model and its generalizations is the dynamical 
mean-field theory (DMFT) <cit.> that has  
provided important insights into multiorbital physics
also in combination with
ab initio calculations for real materials.<cit.>
The continuous-time quantum Monte Carlo method (CTQMC),<cit.> 
particularly the hybridization expansion algorithm (CTHYB),<cit.> is the most widely used impurity solver in multiorbital DMFT calculations. 
However, CTHYB suffers
 from the notorious sign problem when the SOC is included in the calculations.
The sign problem 
grows exponentially with inverse
temperature <cit.> 
and 
typically prevents the study of low-temperature 
symmetry-broken phases.
Alleviating the sign problem in CTHYB would 
help improve our 
understanding of phenomena determined
by the interplay of spin-orbit coupling and correlations.<cit.>

For quantum Monte Carlo (QMC) algorithms based on auxiliary fields, 
there have been various successful advances which
unveiled the origin of the sign problem and suggested a solution in specific
cases,<cit.>
including the recently developed idea of Majorana
symmetry.<cit.>  The rotationally invariant
Hund coupling in the SO-coupled multiorbital Hubbard system,
 however, generates rather complex interaction terms.
Furthermore, the nonlocal-in-time expansion scheme of CTHYB makes it difficult to track the origin of the fermionic sign on the world line configuration. 

In this paper, we systematically study the nature of the sign problem of the
CTHYB for the SO-coupled three-orbital Hubbard model and propose
a strategy to alleviate it.
We employ a numerical sign-optimization scheme, called spontaneous
perturbation stochastic approximation (SPSA),<cit.> to determine the optimal basis 
in terms of average sign. Remarkably, this optimal basis can be well approximated by a simple one,
 which we denote as V_j3/2BA. 
The V_j3/2BA basis is obtained from the relativistic basis by a bonding-antibonding transformation. 



§ MODEL AND METHOD

Our model Hamiltonian is composed of the three terms ℋ = ℋ_t + ℋ_soc + ℋ_int, where ℋ_t, ℋ_soc, and ℋ_int are, respectively,
the electron hopping, SOC, and Coulomb interaction terms.
For the noninteracting electron hopping part we assume a degenerate semicircular density of states ρ(ω) = (1/π D)√(1-(ω/D)^2) and set the half-bandwidth D as the unit of energy.
In the orbital-spin basis V_os, ℋ_soc has the form 

    ℋ_soc = -λ∑^_i∑^_αα'
    σσ'c^†_iασ⟨ασ|𝐋_eff·𝐒|α'σ'⟩ c^_iα'σ' ,

where 𝐋_eff is the l=1 orbital angular momentum operator and 𝐒 is the spin operator.
c^_iασ (c^†_iασ) is the electron annihilation (creation) operator of orbital α (yz,zx,xy) and spin σ (↑,↓) at lattice site i. 
We introduce the Slater-Kanamori form of the Coulomb interaction including the spin-flip and pair-hopping terms: 

    ℋ_int   =    U∑^_i,αn_iα↑n_iα↓ + ∑^_i,α<α' 
    σσ'(U'-J_ Hδ_σσ')n_iασn_iα'σ'
           -J_ H∑^_i,α<α'(c^†_iα↑c^†_iα'↓c^_iα'↑c^_iα↓ + H.c.)
           +J_ H∑^_i,α<α'(c^†_iα↑c^†_iα↓c^_iα'↓c^_iα'↑ + H.c.) .

Here, U (U') is the
intraorbital (interorbital) Coulomb repulsion, and J_H is the strength of the Hund coupling.
To make ℋ_int rotationally invariant, we choose U'=U-2J_H .
This model Hamiltonian covers a wide range of 4d and 5d materials listed in Appendix <ref> by tuning the electron density, λ, U, and J_H .

We solve the model Hamiltonian within the framework of DMFT.<cit.>
The lattice model is mapped onto a quantum impurity model 
by means of a self-consistency relation.
The effective action of the quantum impurity model is written in terms of the local action 𝒮_loc and the hybridization function Δ,

    𝒮_imp   =   𝒮_loc
       +   ∑^_αα'
    σσ'∫_0^βdτ dτ' c^†_ασ(τ)Δ_ασ,α'σ'(τ-τ')c^_α'σ'(τ') .

Here, c^_ασ and c^†_ασ are the fermionic 
annihilation and creation operators 
for orbital α and spin σ for the impurity site.
CTHYB solves this impurity problem by performing an expansion of the partition function in powers of the hybridization function,

    𝒵_imp   =   𝒵_bath∑^∞_k=01/(k!)^2∑^_α_1α'_1
    σ_1σ'_1⋯∑^_α_kα'_k
    σ_kσ'_k
       ×   ∫_0^βdτ_1∫_0^βdτ'_1⋯∫_0^βdτ_k∫_0^βdτ'_k [ Δ^(k)]
       ×   ⟨𝒯_τ c^_α'_1σ'_1(τ'_1)c^†_α_1σ_1(τ_1) ⋯ c^_α'_kσ'_k(τ'_k)c^†_α_kσ_k(τ_k)⟩_𝒮_loc.

Here, Δ^(k) represents the k× k matrix whose elements are 

    Δ^(k)_ij = Δ_α_iσ_i,α'_jσ'_j(τ_i-τ'_j) .

The hybridization function contains the relevant information about the continuous bath degrees of freedom. Since Δ is directly updated in the DMFT self-consistency loop, CTHYB does not rely on a discretization or truncation of the bath, in contrast to other impurity solvers like the exact diagonalization or the numerical renormalization group.<cit.>
In that sense, CTHYB is numerically exact in treating the bath degrees of freedom.

The resulting partition function, Eq. (<ref>), can be expressed as a sum of configuration weights,  

    𝒵_imp = ∑^_𝐱ω(𝐱)

where


    ω(𝐱) = 𝒵_bath(dτ)^2k/(k!)^2⟨𝒯_τ c^_α'_1σ'_1(τ'_1)c^†_α_1σ_1(τ_1) ⋯ c^_α'_kσ'_k(τ'_k)c^†_α_kσ_k(τ_k)⟩_𝒮_loc[ Δ^(k)] .


Here, 𝐱 represents the configuration composed of {k; (α_1,σ_1,τ_1),(α'_1,σ'_1,τ'_1),… (α_k,σ_k,τ_k),(α'_k,σ'_k,τ'_k)}.
One can measure a general time-local observable 𝒪 using MC sampling as

    ⟨𝒪⟩_MC = ∑^_𝐱𝒪(𝐱)ω(𝐱)/∑^_𝐱ω(𝐱) .

A potential sign problem appears when the weight ω(𝐱) becomes negative for certain configurations.
In this case the observable can be sampled through the modified expression 

    ⟨𝒪⟩_MC = ∑^_𝐱sign[ω(𝐱)] 𝒪(𝐱) |ω(𝐱)|/∑^_𝐱sign[ω(𝐱)] |ω(𝐱)| .

A severe sign problem occurs if the small average sign

    ⟨sign⟩=∑^_𝐱sign[ω(𝐱)]|ω(𝐱)|/∑^_𝐱|ω(𝐱)|

amplifies the error propagation in the estimation of the observables.
Especially, for low temperatures, the relative error of the average sign becomes exponentially large as a function of the inverse temperature β.
Thus, reliable estimates of observables cannot be obtained in polynomial time.
An analogous phase problem appears if the configuration weights are complex. 
However, since the average phase is always real, we will for simplicity use the terms sign problem and average sign in the following.


An important aspect to keep in mind is that the sign problem depends on the single-particle basis.
For example, it has been reported that if the electron operators c̃^_j,m are expressed in the relativistic j_eff basis V_jeff,

    (
      [ c̃^_1/2,+1/2; c̃^_3/2,+1/2; c̃^_3/2,-3/2 ])
         =   1/√(6)(
      [  -√(2) +i√(2)  -√(2);     -1      i      2;  -√(3)  i√(3)      0;        ])
      (
      [ c^_yz,↓; c^_zx,↓; c^_xy,↑;         ]) ,
    (
      [ c̃^_1/2,-1/2; c̃^_3/2,-1/2; c̃^_3/2,+3/2 ])
         =   1/√(6)(
      [  -√(2) -i√(2)  +√(2);     +1      i      2;  +√(3)  i√(3)      0;        ])
      (
      [ c^_yz,↑; c^_zx,↑; c^_xy,↓;         ]),

the average sign is improved <cit.> due to the diagonalized hybridization function matrix.<cit.>  (An average sign of 1 means no sign problem, while an average sign approaching 0 means a severe sign problem.) 

To further optimize the basis beyond V_jeff, we introduce the
numerical optimization scheme SPSA. In this scheme the average sign becomes the
objective function on the parameter space where each point corresponds to a
single-particle basis used in CTHYB.
The SPSA approximates the gradient for a given parameter point.
At each iteration, the objective function is measured at both the positively and negatively perturbed parameter points along a stochastically chosen direction.
The gradient is evaluated from this two-point measurement, and the parameter point is updated.
Compared to the finite-difference stochastic approximation which involves a number of measurements proportional to the dimension of the parameter space, the SPSA potentially converges faster in time if the parameter space is high-dimensional and the measurement of the objective function is computationally expensive.





§ RESULTS

In Fig. <ref> we plot the average sign for the four different bases considered in this work as a function of electron density
[Fig. <ref>(a)]  and spin-orbit coupling strength [Fig. <ref>(b) and (c)]. 
The orbital-spin basis V_os has a severe sign problem 
since a drastic drop in the average sign of V_os appears as a function of the SOC strength both below and above half-filling. 
Figure <ref> shows the evolution of the average sign as a function of U (J_H) at various electron fillings.
The suppressed average sign of V_os in the noninteracting limit 
  implies that the source of the sign problem in this basis is the off-diagonal hybridization function generated by the SOC.
The basis V_jeff, defined in Eq. (<ref>), diagonalizes
the hybridization-function matrix of the Hamiltonian and thus recovers the
sign-free behavior in the noninteracting limit, as demonstrated in Fig. <ref>.  In the presence of a nonzero
Hund coupling J_H, however, V_jeff shows a sign problem
as well, especially at intermediate interaction strengths and away from half-filling.





In order to further improve the average sign we introduce the basis V_j3/2BA in which c̃^_3/2,+3/2 and c̃^_3/2,-3/2 are mixed in the bonding-antibonding manner:

    (
      [ ã^_3/2,+; ã^_3/2,- ])
      = 1/√(2)(
      [  1  1;  1 -1;    ])
      (
      [ c̃^_3/2,+3/2; c̃^_3/2,-3/2 ]) .

A similar transformation was previously explored for an isolated trimer.<cit.> It turns out that V_j3/2BA is superior to V_jeff for most parameter sets investigated.
In the noninteracting limit, V_j3/2BA (like V_jeff) yields an average sign of unity (see Fig. <ref>). 
As U is switched on, V_j3/2BA results in a larger average sign than V_jeff .
In Fig. <ref>, the improvement of the average sign in V_j3/2BA is
 prominent in the itinerant regime with intermediate Coulomb interaction strength.
This improvement persists over the whole range of electron densities for U=2 as shown in Fig <ref>(a) .
In the Mott insulating regime (marked by open squares and circles in Figs. <ref> and <ref>), 
on the other hand, the difference between the average signs of
V_j3/2BA and V_jeff is smaller.  Since the CTHYB
is based on the expansion around the localized limit, those observations imply
that V_j3/2BA effectively prevents sign-problematic high-order
processes.

Note that neither V_j3/2BA nor V_jeff is an optimal basis for the non-SO-coupled model.
As we show in Figs. <ref>(b) and (c), at the
non-SO-coupled point (λ=0), V_j3/2BA and V_jeff
exhibit a sign problem, in contrast to V_os .
This demonstrates that Eq. (<ref>) includes dangerous interacting terms in the V_j3/2BA and V_jeff basis.
One of these terms is the correlated-hopping (CH) term, which has the form 

    √(2)J_H/3∑^_s=±[s(2ñ_3/2,s3/2 - ñ_3/2,s̅3/2- ñ_3/2,s̅1/2)c̃^†_1/2,s1/2c̃^_3/2,s1/2 + H.c.]

in V_jeff and 

    √(2)J_H/6∑^_s=±[ s(ñ_3/2,+ + ñ_3/2,-   - 2 ñ_3/2,s̅1/2) c̃^†_1/2,s1/2c̃^_3/2,s1/2+H.c.]
          +J_H/2[(ñ_1/2,+1/2-ñ_1/2,-1/2)ã^†_3/2,+a^_3/2,- + H.c.]

in V_j3/2BA .
This CH term is the major source of the sign problem in the SO-coupled Hamiltonian.
Table <ref> analyzes the effect of different terms in the local
 Hamiltonian on the average sign for both V_j3/2BA and V_jeff .
There is a substantial drop in the average sign when the CH terms are introduced.





In V_j3/2BA and V_jeff, there emerge other new terms involving four different flavors, which are not of the spin-flip type or pair-hopping type appearing in the Slater-Kanamori Hamiltonian in V_os .
Those terms are denoted as four-scattering terms in Table <ref> .
The average sign in both V_j3/2BA and V_jeff becomes even lower when the CH term is combined with the pair-hopping and the four-scattering terms.
The full interacting Hamiltonian in V_jeff is explicitly written in Appendix <ref>.
Furthermore, we analyze the nature of the self-consistent solutions for the Hamiltonians without the CH or FS terms and we discuss the potential use of such masked Hamiltonians in Appendix <ref>.




In what follows we will determine by the SPSA method the optimal basis in terms of the average sign for the
parameters with the most severe sign problem and show
that this optimal basis is nearly identical to the V_j3/2BA basis.
For that,
we search the basis space generated by the SO(4) rotation group for j=3/2, whose 4× 4 matrix representation is denoted by M and transforms the electron operators in V_jeff as follows: 

    ã^_3/2 = M·c̃^_3/2 ,
    ã^⊺_3/2 = (ã^_3/2,1,ã^_3/2,2,ã^_3/2,3,ã^_3/2,4) ,
    c̃^⊺_3/2 = (c̃^_3/2,+1/2,c̃^_3/2,-1/2,c̃^_3/2,+3/2,c̃^_3/2,-3/2) .

Since the off-diagonal hybridization function is a clear source of the severe sign problem as shown in the cases of V_os (Figs. <ref> and <ref>), we exclude the mixing between the j=1/2 and 3/2 subspaces to preserve the diagonal structure of the hybridization.
Without mixing between the j=1/2 and 3/2 subspaces, one can fix the basis for the j=1/2 subspace using the rotational symmetry generated by the total angular momentum operator, 𝐉_eff=𝐋_eff+𝐒 without loss of generality.


To parametrize the basis space, we introduce the isoclinic decomposition of M as
	M = M_LM_R ,
where

    M_L=(
    	[  a -b -c -d;  b  a -d  c;  c  d  a -b;  d -c  b  a ]),M_R=(
    	[  p -q -r -s;  q  p  s -r;  r -s  p  q;  s  r -q  p ]).

Here, a^2+b^2+c^2+d^2 = 1 and p^2+q^2+r^2+s^2 = 1 .
Under these constraints among {a,b,c,d} and {p,q,r,s} the dimension of the parameter space becomes 6.

Figure <ref> shows how the SPSA works while searching for the optimal basis in this parameter space.
The evolution of the average sign as a function of number of iterations
at both positively (⟨sign^+⟩) and negatively (⟨sign^-⟩) perturbed points in the parameter space is shown in Fig. <ref> (a).
The average sign value converges to ∼ 0.74 for n=3.5, T=0.03, λ=0.25, U=2.0, and J_H/U=0.15 .
Within numerical accuracy, it is very close to the value of
V_j3/2BA defined in Eq. (<ref>) .  This shows that the
V_j3/2BA basis is at least near the local optimum in parameter
space. 
The inset of Fig. <ref>(a) illustrates the convergence of the parameters in Eqs. (<ref>). 

Figure <ref>(b) and (c) show the SPSA sequence in a small parameter subspace. 
We introduce two parameters, θ_1 and θ_2, representing the restricted basis transformation 

    M
    	= (
    	[  cosθ_1  sinθ_1       0       0; -sinθ_1  cosθ_1       0       0;       0       0  cosθ_2  sinθ_2;       0       0 -sinθ_2  cosθ_2;         ]) .

 Figure <ref>(c) shows that the landscape of the average sign is smooth, so that the
SPSA search based on the gradient approximation can successfully find the local optimum in that subspace.
Furthermore, V_j3/2BA, corresponding to θ_1=0 and θ_2=π/4, is shown to be very close to the optimum found by the SPSA.
Figure <ref>(b) plots the trajectory determined by the SPSA.



We finally investigate the temperature scaling of the average sign for the different bases.
Figure <ref> shows the exponential scaling ⟨sign⟩∼exp(-βδ F + F_0) 
as a function of the inverse temperature β. 
Here, δ F can be regarded as the free-energy difference from the auxiliary bosonic system without a sign problem, which determines the temperature scaling. 
For both fillings, n=3.5 and 2.5, the V_j3/2BA basis shows an improved temperature scaling exponent δ F and a comparable offset F_0 .



§ CONCLUSION

We have investigated the nature of the sign problem in the CTHYB for the SO-coupled multiorbital Hubbard model.
We found that the correlated hopping term that
appears in the j_eff basis is the major source of the sign problem,
 and we introduced a basis–the j=3/2 bonding-antibonding basis–which alleviates the effects of those terms.
By applying the stochastic optimization scheme, we found that
the j=3/2 bonding-antibonding basis is near a local optimum in the parameter space considered.
These results (i) provide useful guidelines for the choice
 of the single-particle basis in CTHYB simulations of SO-coupled systems,
(ii) introduce an algorithm to numerically determine the optimal basis, and 
(iii) suggest using the intrinsic fermionic sign bound in a
systematic search of the full SO(4) basis space
as a strategy for the further improvement of CTHYB for this class of systems. 

In addition, our findings have relevance for descriptions beyond DMFT that include
nonlocal correlations. For example, 
since the average sign is alleviated at the level of the partition function in V_j3/2BA,
we expect our results to be helpful for a reliable estimation of two-particle correlation functions,
which are the essential building blocks of diagrammatic extensions of DMFT.<cit.> 
For cluster extensions of DMFT, on the other hand, our numerical basis optimization scheme can be applied to an extended basis space including site (cluster) degrees of freedom. 
A promising strategy might be to combine the V_j3/2BA basis with a bonding-antibonding transformation that breaks loops on the cluster.<cit.>

Note added. During the final stages of this work, a series of conceptually related studies appeared on  arXiv.<cit.>


A.J.K. was supported by EPSRC through Grant No. EP/P003052/1. 
P.W. acknowledges support from the Swiss National Science Foundation via NCCR Marvel. 
R.V. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) 
through Grant No. VA117/15-1. 
The computations were performed at the Centre for Scientific Computing (CSC)
in Frankfurt.







§ ENERGY SCALES IN MATERIALS


In this appendix, we list materials related to our spin-orbit-coupled (SOC) three-orbital model, and the corresponding energy scales.  
The SOC is an intrinsic relativistic effect of the atom whose strength is an increasing function of the atomic number.
Hence, 5d materials (Ta, W, Os, and Ir) consistently show a larger SOC strength compared to 4d materials (Ru and Rh).
Table <ref> lists materials with various electron configurations which correspond to the electron density in our model.
The ratio between the SOC strength and the half-bandwidth used in our model is given in the fourth column.


Since the proposed V_j3/2BA basis alleviates the sign problem in a wide parameter range (electron density and strength of the SOC and many-body interaction), as shown in Figs. 1 and 2 of the main text, it should enable improved material-specific calculations for a broad range of materials.




§ INTERACTING HAMILTONIAN IN THE V_JEFF BASIS


In Eq. (<ref>), we classify the interacting Hamiltonian in the V_jeff basis into five classes: density-density (ℋ_dd), correlated-hopping (ℋ_ch), pair-hopping (ℋ_ph), four-scattering (ℋ_fs), and spin-flip (ℋ_sf) terms.
The Hamiltonian in the V_j3/2BA basis is classified in the same way.


    ℋ_int   =   ℋ_dd + ℋ_ch + ℋ_ph + ℋ_fs + ℋ_sf ,
    ℋ_dd   =    (U- 4J_H/3)ñ_k/2,+1/2ñ_1/2,-1/2 + (U-J_H)(ñ_3/2,1/2ñ_3/2,-1/2 + ñ_3/2,3/2ñ_3/2,-3/2) 
          + (U-7J_H/3)∑^_s=±(ñ_1/2,s1/2ñ_3/2,s̅1/2 + ñ_3/2,s1/2ñ_3/2,s3/2 + ñ_3/2,s1/2ñ_3/2,s̅3/2)
          + (U-2J_H)∑^_s=±ñ_1/2,s1/2ñ_3/2,s1/2 + (U-8J_H/3)∑^_s=±ñ_1/2,s1/2ñ_3/2,s̅3/2 + (U-5J_H/3)∑^_s=±ñ_1/2,s1/2ñ_3/2,s3/2 ,
    ℋ_ch    =    √(2)J_H/3∑^_s=±[s(2ñ_3/2,s3/2 - ñ_3/2,s̅3/2- ñ_3/2,s̅1/2)c̃^†_1/2,s1/2c̃^_3/2,s1/2 + H.c.] ,
    ℋ_ph   =   - 5J_H/3[c̃^†_1/2,+1/2c̃^†_1/2,-1/2(c̃^_3/2,-1/2c̃^_3/2,+1/2 - c̃^_3/2,-3/2c̃^_3/2,+3/2) + H.c.] - 4J_H/3[c̃^†_3/2,+1/2c̃^†_3/2,-1/2c̃^_3/2,-3/2c̃^_3/2,+3/2 + H.c.] ,
    ℋ_fs   =   - J_H/√(3)∑^_s=±[c̃^†_1/2,s1/2c̃^†_3/2,s1/2c̃^_1/2,s̅1/2c̃^_3/2,s3/2  + c̃^†_1/2,s1/2c̃^†_3/2,s̅3/2c̃^_1/2,s̅1/2c̃^_3/2,s̅1/2 + H.c.]
          - √(2)J_H/√(3)∑^_s=±s[c̃^†_1/2,s1/2c̃^†_3/2,s̅1/2c̃^_3/2,s̅3/2c̃^_3/2,s3/2 + c̃^†_1/2,s1/2c̃^†_3/2,s1/2c̃^_3/2,s̅1/2c̃^_3/2,s3/2 + H.c.] ,
    ℋ_sf   =   2J_H/3∑^_s=±[c̃^†_1/2,s1/2c̃^†_3/2,s̅1/2c̃^_3/2,s1/2c̃^_1/2,s̅1/2 + H.c.] .





§ SELF-CONSISTENT SOLUTION OF THE MASKED HAMILTONIANS


In this appendix, we investigate the nature of the self-consistent solutions for masked Hamiltonians.
Since we can increase the average sign substantially by dropping the most problematic correlated-hopping (CH) or four-scattering (FS) terms, 
those masked Hamiltonians potentially provide a useful approximation of the full Hamiltonian if the self-consistent solution is sufficiently close to the one for the full Hamiltonian.

The masking of those terms, however, modifies the local eigenstates.
Table <ref> shows the form of the ground states of the full local Hamiltonian in V_jeff defined in the main text.
When we mask the CH or FS terms, the ground-state degeneracy of the full local Hamiltonian is broken in the N=2 and 3 sectors.
Tables <ref> and <ref> present the ground states for the masked local Hamiltonians without the CH and FS terms, respectively.
The remaining ground states for the N=2 and 3 sectors depend on the type of the masked terms.
When the FS terms are dropped, the form of the highest |J_z| ground states remains the same as that for the full Hamiltonian.
Since the FS terms involve four different flavors by definition, they become irrelevant for the highest |J_z| states with fixed j=3/2 and m_j=± 3/2 electrons.
On the other hand, masking the CH terms selects the |J_z|=1 states for the N=2 sector and the |J_z|=1/2 states for the N=3 sector with slightly modified coefficients, which demonstrates the relevance of these terms for the highest |J_z| states.

Such a modified degeneracy of the local Hamiltonian leads to substantial changes in observables, especially when the system becomes localized.
In Figs. <ref>(a), <ref>(c), and <ref>(e), for example, one can see a sizable difference between the electron densities from the full and masked local Hamiltonians for the large-U Mott insulator.
The degeneracy of the electron density between the j=3/2, m_j=± 1/2, and m_j=± 3/2 flavors is naturally broken for the masked Hamiltonian.
Mott localization is signaled by the suppression of the spectral function at the Fermi level, and this quantity can be approximately evaluated as Ã(ω=0)=-1/π T G(τ=β/2) .
Moreover, the Mott transition point U_c for n=2 and 4 is reduced as we mask the CH or FS terms.
Compared to the FS-dropped Hamiltonian, the CH-dropped one shows a further reduction in U_c for the n=2 case.
This kind of U_c reduction as a result of a degeneracy breaking was reported in a non-spin-orbit-coupled system in the presence of a single-particle crystal-field splitting <cit.> and in the absence of the spin-flip and pair-hopping terms (so-called Ising-type Hund coupling) at the many-body level <cit.> (see Fig. <ref>).


One interesting observation is that the density values of the full Hamiltonian are approximately recovered by the masked one if we artificially symmetrize the j=3/2, m_j=± 1/2, and m_j=± 3/2 flavors during the self-consistency loop of the DMFT.
Figure <ref> shows the corresponding electron density and the approximated spectral function at the Fermi level. The density values become much closer to those of the full Hamiltonian.
Especially, the FS-dropped Hamiltonian at and below half-filling yields results which are consistent within ∼ 5% relative error.
As a result of the modification of the Hamiltonian, the sign problem of the CTHYB simulation is alleviated by up to ∼ 15%. When applied with proper care, this symmetrization trick could provide a useful estimate for physical observables when the original Hamiltonian cannot be treated due to the serious sign problem.















