Tackling the Sign Problem in the Doped Hubbard Model with Normalizing Flows

T ac kling the Sign Problem in the Dop ed Hubbard Mo del with Normalizing Flo ws Dominic Sc huh , 1, 2 , ∗ Lena F unc ke , 1, 2 Janik Kreit , 1, 2 Thomas Luu , 2, 3 and Simran Singh 1, 2 1 T r ansdisciplinary R ese ar ch Ar ea (TRA) Matter, University of Bonn, Bonn, Germany 2 Helmholtz Institute for R adiation and Nucle ar Physics (HISKP), University of Bonn, Bonn, Germany 3 Institute for A dvanc e d Simulation 4 (IAS-4), F orschungszentrum J¨ ulich, J¨ ulich, Germany (Dated: March 20, 2026) The Hubbard mo del at finite chemical potential is a cornerstone for understanding dop ed correlated systems, but sim ulations are sev erely limited b y the sign problem. In the auxiliary-field form ulation, the spin basis mitigates the sign problem, y et sev ere ergo dicity issues ha ve limited its use. W e extend recen t adv ances with normalizing flows at half-filling to finite c hemical potential by introducing an annealing scheme enabling ergo dic sampling. Compared to state-of-the-art h ybrid Monte Carlo in the c harge basis, our approach accurately repro duces exact diagonalization results while reducing statistical uncertain ties by an order of magnitude, opening a new path for simulations of dop ed correlated systems. Intr o duction —Understanding systems of strongly cor- related electrons is fundamental to quantum man y-bo dy- ph ysics [1]. One of the most widely used theoretical framew orks for studying suc h systems is the Hubbard mo del [2 – 4], whic h captures the essential comp etition b et w een electron kinetic energy and on-site in teractions, thereb y giving rise to phenomena suc h as the Mott in- sulating b ehavior [5], magnetism and sup erconductiv- it y [6 – 8]. Over the years, a wide range of metho ds ha ve b een developed to analyze the Hubbard mo del. In the w eak-interaction regime, perturbative approaches pro- vide v aluable insigh ts [9]. In contrast, in the strongly correlated regime, where p erturbative metho ds fail, non- p erturbativ e metho ds are required. These can be broadly classified into Hamiltonian-based approac hes and action- based approaches. The latter typically rely on the finite- temp erature auxiliary-field form ulation of the Hubbard mo del, where Monte Carlo simulations b ecome essential (see, e.g., Refs. [10 – 17] and references therein). Dop ed systems are describ ed by a finite chemical p o- ten tial and are essential for describing realistic materials and molecular systems. How ev er, a wa y from half-filling or on non-bipartite lattices, the Hubbard mo del suffers from a severe numerical sign problem, making reliable es- timates of observ ables increasingly challenging [18]. F ur- thermore, due to practical ergo dicity problems in the spin basis, Monte Carlo based metho ds are usually restricted to work in the charge basis of the Hubbard mo del, where the sign problem is particularly pronounced. Although sev eral metho ds ha v e b een developed to mitigate these c hallenges [19, 20], they still fall short for crucial regions of the parameter space. In this work, we build on the results of Refs. [21, 22] and dev elop the first deep generativ e mac hine learning approac h to the finite-temperature auxiliary-field formu- lation of the Hubbard mo del at non-zer o c hemical p oten- tial. T o this end, w e in troduce an annealing sc heme to ov er- come the practical ergo dicity problem in the spin basis [23] and benchmark our metho d against results obtained b y state-of-the-art optimized h ybrid Monte Carlo (HMC) sim ulations [20]. The remainder of this letter is organized as follo ws. W e first in tro duce the H ubbard mo del in the spin basis and highligh t its k ey features, follow ed by a brief ov erview of normalizing flo ws. W e then present our prop osed anneal- ing scheme and its results, b efore concluding. The Hubb ar d mo del at finite density —The Hubbard mo del [2 – 4, 24] provides a simple framework for study- ing the interaction of electrons in materials and thereby their conducting prop erties. In its simplest form, the dy- namics consist of nearest-neighbor hopping, a repulsive Coloum b-like on-site in teraction, and a chemical p oten- tial. In this work, we will fo cus on studying the mo del in the spin basis and at finite chemical p oten tial. The corresp onding Hamiltonian can b e written as H = − κ X ⟨ x,y ⟩ ,σ c † x,σ c y ,σ − U 2 X x ( n x, ↑ − n x, ↓ ) 2 − µ X x q x , (1) where κ is the hopping parameter, U the on-site in ter- action, and µ the c hemical p oten tial. Here, c † x,σ ( c x,σ ) denotes the creation (annihilation) op erator for an elec- tron at site x with spin σ (= ↑↓ ), n x,σ = c † x,σ c x,σ is the corresp onding n um b er op erator, and q x = P σ n x,σ is the net charge at site x . T o circumv en t the exponential gro wth of the Hilb ert space with the n umber of spatial sites N x , the mo del can b e reformulated in an action-based auxiliary field represen tation [2 – 4, 24–26], enabling the use of Monte- Carlo based metho ds. In this framework, observ ables are expressed as statistical quantum-mec hanical expectation v alues, ⟨O ⟩ = 1 Z tr  e − β H O  , (2) where the partition function is Z = tr  e − β H  , the in- v erse temp erature is β = 1 /k B T , and the trace runs o ver the entire F o ck space. 2 After applying a Suzuki-T rotter decomposition [27, 28] and a Hubbard-Stratono vic h transformation [29] to de- couple the quartic fermion in teraction, this exp ectation v alue can b e rewritten in terms of an auxiliary b osonic field ϕ , ⟨O ⟩ = 1 Z Z D [ ϕ ] e − S [ ϕ ] O [ ϕ ] . (3) Here, the Hubbard action S [ ϕ ] is given by S [ ϕ ] = 1 2 ˜ U X x,t ∈ Λ ϕ 2 x,t − ln det ( M [ ϕ, ˜ κ, ˜ µ ] M [ − ϕ, ˜ κ, ˜ µ ]) , (4) where Λ = N x × N t is the lattice v olume, ϕ xt denotes the field at site ( x, t ), and ˜ U = U β / N t , ˜ κ = κβ / N t , ˜ µ = µβ / N t are the dimensionless lattice couplings. The fermion matrices M [ ± ϕ, ˜ κ, ˜ µ ] originate from the spin-up and spin-down electrons in Eq. (1), resp ectively . In this reform ulation, one can choose a discretization of the fermion matrix [17]. In this work, we adopt the exp onen tial discretization, as it preserves c hiral (or sub- lattice) symmetry [23, 30]. In this discretization, the ma- trix elements are given by M [ ϕ | ˜ κ, ˜ µ ] x ′ t ′ ,xt = δ x ′ ,x δ t ′ t − [ e h ] x ′ x e ϕ xt + ˜ µ B t ′ δ t ′ ,t +1 . (5) Here, h is the hopping matrix defined as h = ˜ κδ ⟨ x,x ′ ⟩ , with ⟨ x, x ′ ⟩ denoting nearest-neighbor-interactions, and B incorp orates anti-perio dic b oundary conditions in the temp oral direction [31]. The observ able considered in this work is the single- particle Euclidean-time correlator, defined as C x,y ( t ) ≡ ⟨ c x ( t ) c † y (0) ⟩ = 1 N cfg X { ϕ } M − 1 ( t,x );(0 ,y ) [ ϕ | ˜ κ, ˜ µ ] . (6) These correlators are computed for all site pairs ( x, y ) and diagonalized to obtain the eigenfunctions C ± ( t ) [22, 23], sampled on N cfg configurations. This observ able has previously been used to illustrate the lac k of ergo dict y in sampling the Boltzmann distribution [22, 23]. In this w ork, w e demonstrate the ergodic nature of our no v el annealing-based sampling method by comparing these correlators with results obtained from exact diagonaliza- tion (ED), where av ailable. Choic e of b asis —The Hubbard mo del is most com- monly sim ulated either in the repulsiv e spin basis [15, 32] or the attractive charge basis [17, 33]. In the action-based form ulation, the difference betw een the t w o lies in how the auxilliary field ϕ xt en ters the fermion matrix in Eq. 5. In the spin basis the fields are real-v alued, making the matrix elements real, whereas in the charge basis, the fields are purely imaginary . Mon te Carlo simulations of the Hubbard mo del typi- cally face tw o ma jor challenges: (i) ergo dicit y issues [15, 23] and (ii) the non-p ositivit y of the probabilit y weigh t arising from the pro duct of fermion determinan ts [34, 35]. Consequen tly , at ze ro chemical p otential, the feasibilit y of HMC sim ulations dep ends primarily on effectively ad- dressing the ergo dicity problem. This aspect has been systematically in vestigated in Refs. [23, 36], where the authors sho w that the formal ergo dicit y issues of the c harge basis are often easier to handle than the practical ergo dicit y issues of the spin basis. Since the primary ob jective of this w ork is to inv esti- gate the system at finite chemical p oten tial, it is adv an- tageous to adopt a formulation that offers b etter con trol o ver the sign problem. In this context, as we explicitly demonstrate, the spin basis offers a structural adv an tage: the fermion matrix is real for each spin species separately , thereb y a voiding the complex probabilit y w eigh ts that arise in the charge basis. As a result, the sign problem is purely real rather than complex, which simplifies its n umerical treatment. In the spin basis, the fermion matrices of the tw o sp ecies are real but not related b y conjugation, so the asso ciated probabilit y weigh t is not manifestly p ositiv e. This can lead to a p otential sign problem ev en at half- filling. How ev er, as shown in Ref. [23], the fermion deter- minan t in the exp onen tial discretization satisfies a sym- metry relation b et ween the tw o sp ecies that guarantees p ositivit y of the probability w eight at zero chemical p o- ten tial. In tro ducing a finite chemical p otential explicitly breaks this symmetry and thus destroys the exact p osi- tivit y of the measure [37]. Sign quenching & r e-weighting— When sim ulating the- ories with a complex action, standard sampling-based metho ds — whic h rely on a probabilistic interpretation of the Boltzmann w eight e − β S — are no longer directly applicable. In such cases, one typically follows one of tw o strategies: (i) simulate the phase-quenche d theory [38], in which the complex phase of the Boltzmann weigh t is neglected and subsequen tly restored by rew eighting; or (ii) simulate the sign-quenche d theory [39, 40], where the Boltzmann weigh t is replaced b y the absolute v alue of its real part, and sign fluctuations are rein tro duced through sign r eweighting . The phase-quenc hed approach is pri- marily hindered b y severe o v erlap problems [41], whereas the sign-quenched approac h suffers from an exp onentially degrading signal-to-noise ratio [42, 43]. In the spin basis, the sign problem is liter al ly a sign problem: the probability weigh t is real but not p ositive definite, rather than complex. W e therefore emplo y sign quenc hing, defining w q = | w t | for real weigh ts w t that fluctuate in sign. Observ ables in the full theory , ⟨ O ⟩ t , are recov ered via reweigh ting: ⟨ O ⟩ t = R D ϕ w t ( ϕ ) O ( ϕ ) R D ϕ w t ( ϕ ) = ⟨ w t w q · O ⟩ q ⟨ w t w q ⟩ q , (7) where ⟨·⟩ q denotes sampling from the sign-quenc hed the- ory . 3 The sev erity of the sign problem is quantified by the denominator in Eq. (7), commonly referred to as the av- er age sign Σ. A v alue of Σ = 1 corresp onds to the ab- sence of a sign problem, whereas Σ → 0 indicates the most sev ere sign problem, where cancellations b etw een configurations suppress the signal. The resulting reduc- tion in statistical efficiency leads to an effectiv e sample size N eff = | Σ | 2 N [44], with Monte Carlo uncertainties scaling as ( N eff ) − 1 / 2 . With this in mind, w e mo del the sign-quenche d distri- bution e − β S q , where S q [ ϕ ] = X x,t ∈ Λ ϕ 2 x,t 2 ˜ U − log | det ( M [ ϕ | ˜ κ, ˜ µ ] M [ − ϕ | ˜ κ, ˜ µ ]) | . (8) This gives us a w ell-defined probabilit y weigh t, from whic h the physical theory can b e recov ered by rew eight- ing observ ables with the sign of the original weigh t. Normalizing flows — Normalizing flows (NFs) [45, 46] are a class of generativ e mo dels that enable efficien t sampling from high-dimensional distributions. This is ac hieved b y learning a bijective map f θ b et w een a simple prior distribution q z and a target distribution q y , f θ : z ∼ q z 7→ y ∼ q y , (9) where f θ is parametrized by neural netw ork w eigh ts θ . In our setup, the target degrees of freedom y are the aux- iliary field configurations of the Hubbard mo del. Con- structing these maps requires balancing the expressivit y of the neural net work with the need for a tractable Ja- cobian. Several approaches hav e b een developed to this end, such as coupling-based NFs [47 – 49], autoregressiv e NFs [50], and contin uous NFs [51]. In this work, we employ a coupling-based NF with a real-v alued non-volume preserving (RealNVP) archi- tecture [48]. Sp ecifically , we consider a sequence of L smo oth, in vertible mappings f l θ l : R | Λ | → R | Λ | , each de- fined by a neural net work with parameters θ l . Since each la yer is bijective, their comp osition defines an ov erall in- v ertible transformation: y ≡ f θ ( z ) =  f L θ L ◦ f L − 1 θ L − 1 ◦ · · · ◦ f 1 θ 1  ( z ) , (10) where θ = { θ 1 , . . . , θ L } denotes the full set of train- able parameters. A crucial prop erty of this construc- tion is that the determinan t of the Jacobian of f θ can b e computed efficiently , which is required for tractable lik eliho o d-based training and density estimation of q y . The generated field configurations ϕ are then given by ϕ ≡ f θ ( z ) , (11) with the asso ciated learned density q θ ( f θ ( z )) = q z ( z ) ·     ∂ f θ ∂ z     − 1 . (12) T raining of the neural net work pro ceeds by minimiz- ing the rev erse Kullbac k–Leibler (KL) div ergence [52], whic h measures the distance b etw een the flow-induced distribution q θ and the target Boltzmann distribution p ( ϕ ) = e − S [ ϕ ] / Z . This ob jective can b e written as KL( q θ || p ) = E z ∼ q z h log q θ ( f θ ( z )) + S [ f θ ( z )] i , (13) where S [ · ] denotes the action of the underlying phys- ical system. T o ensure exact and unbiased sampling, observ ables calculated from the flow are subsequen tly rew eighted according to the unnormalized imp ortance w eights w ( ϕ ) = exp( − S [ ϕ ]) /q θ ( ϕ ), where, in our case, S = S q [ ϕ ] corresp onds to the sign-quenched action. Avoi ding pr actic al er go dicty pr oblems with normaliz- ing flows —It is w ell kno wn that training with the re- v erse KL divergence can lead to mo de dropping in the learned densit y [53 – 55]. In challenging settings, suc h as the Hubbard mo del in the spin basis, having an ergo dic sampler is essen tial to ensure un biased results [22, 23]. Symmetry-informed architectures, e.g. canon- icalization [21] or symmetry-enforcing sto chastic mo du- lation [22, 56], hav e prov en highly effective at sampling distributions with a strong practical ergo dicity problem. Ho wev er, their applicability is limited to systems exhibit- ing particular symmetries, which in the Hubbard mo del t ypically arise only in the strongly coupled regime, where these symmetries b ecome particularly pronounced. T o describ e systems across a range of couplings, w e cannot rely on deterministic symmetry transformations. Instead, w e emplo y an annealing scheme [57], in which an auxiliary parameter λ scales the fermionic contribution to the action, S q ,λ [ ϕ ] = X x,t ∈ Λ ϕ 2 x,t 2 ˜ U − λ · log | det ( M [ ϕ | ˜ κ, ˜ µ ] M [ − ϕ | ˜ κ, ˜ µ ]) | . (14) F or λ = 0, the fermion determinant is absent, leaving a simple Gaussian distribution, while λ = 1 recov ers the full interacting Hubbard action. Gradually increasing λ during training smo othly deforms the target distribu- tion tow ard the final multimodal target, as illustrated sc hematically in Fig. 1. This mitigates the mo de-seeking b eha vior of the reverse KL divergence, allowing the NF to receiv e training signals from all mo des and represent ev en widely separated probability sectors. This approach is highly scalable, as its computational cost grows only moderately with system size, adding a minimal training o verhead. W e emphasize, how ev er, that while the annealing pro cedure emprically demonstrates strong robustness against mode dropping, it do es not pro vide a formal guarantee of complete mo de cov erage in the final mo del. R esults —In the follo wing, we ev aluate the performance of our annealing-based NF ansatz for the Hubbard model in the spin basis. As a benchmark, we use state-of-the 4 λ = 0 λ = 0 . 5 λ = 1 FIG. 1. Schematic illustration of the annealing pro cedure controlled b y the parameter λ . At λ = 0, the fermion determinant is absent and the target distribution is Gaussian. As λ is increased, fermionic effects are gradually introduced, deforming the distribution until the full interacting Hubbard action with its multimodal structure is recov ered at λ = 1. art HMC in the charge basis with optimized shifts [20]. W e compare the av erage sign of both approac hes and calculate correlation functions for the most c hallenging system. Each model is trained for 50k steps, with the an- nealing parameter λ increased linearly from 0 to 1 during the first 1k steps. Figure 2 shows the av erage sign obtained with the NF ansatz compared to optimized HMC simulations for an eigh t-site hexagonal lattice with U = 2, β = 8, and N t = 16. Although optimized shifts are known to al- leviate the sign problem in the c harge basis compared to standard HMC [20], w e still observe a significan t suppres- sion of the av erage sign for µ ∈ [1 . 0 , 2 . 0] in this setup. In con trast, simulations p erformed in the spin basis using the NF approac h main tain a substantially larger a verage sign across this region, with an improv emen t of approxi- mately a factor of four at the minimum. 0 1 2 3 4 5 6 C h e m i c a l P o t e n t i a l 0.2 0.4 0.6 0.8 1.0 A v e r a g e S i g n O+HMC CH NF SP FIG. 2. Average sign Σ as a function of the c hemical potential µ for an eight-site hexagonal lattice with U = 2, β = 8, and N t = 16. Sho wn are results from the annealing-based NF ap- proac h in the spin basis (NF SP , yello w) and optimized HMC sim ulations in the charge basis [20] (O+HMC CH, grey). T o v alidate and assess our approach, we calculate one- b ody correlation functions C ( t ), defined in Eq. 6, and compare them to the ground truth obtained via ED as w ell as the state-of-the-art optimized HMC algorithm. W e consider an eight-site hexagonal lattice with param- eters U = 2, β = 8, µ = 1 . 75, and N t = 40, whic h represen ts the largest exactly solv able system exhibiting the most severe sign problem, as indicated in Fig. 2. The corresp onding results are shown in Fig. 3. W e find that the NF approach repro duces the exact results with high precision while reducing uncertainties by an order of mag- nitude compared to HMC. Additional results for U = 3 are provided in App. C. Lastly , we extend our analysis beyond exactly solv- able systems and compute one-b o dy correlation func- tions for an 18-site hexagonal lattice with parameters 0 2 4 6 8 t 1 0 4 1 0 3 1 0 2 1 0 1 1 0 0 C ( t ) ED O+HMC CH NF SP FIG. 3. One-b ody correlation functions for an eight-site hexagonal lattice with U = 2, β = 8, µ = 1 . 75, and N t = 40. Sho wn are continous exact diagonalization results (ED, grey dashed), state-of-the-art optimized HMC in the charge basis (O+HMC CH, grey), and our annealing-based NF ansatz in the spin basis (NF SP , yello w). Both HMC and NF results are generated with 1M samples eac h. 5 0 2 4 6 8 t 1 0 5 1 0 4 1 0 3 1 0 2 1 0 1 1 0 0 C ( t ) O+HMC CH NF SP FIG. 4. One-bo dy correlation functions for an 18-site hexag- onal lattice with U = 2, β = 8, µ = 2 . 5, and N t = 40. Sho wn are state-of-the-art optimized HMC in the charge ba- sis (O+HMC CH, grey) and our annealing-based NF ansatz in the spin basis (NF SP , yello w). Both HMC and NF results are generated with 1M samples eac h. U = 2 , β = 8 , µ = 2 . 5 , and N t = 40, as sho wn in Fig. 4. In the absence of exact results at this system size, w e v alidate our approac h using a less noise-sensitive corre- lator, finding go o d agreemen t with HMC (top correlator in Fig. 4). F or more noise-sensitive observ ables (b ottom correlator), the NF metho d pro duces smaller uncertain- ties and impro ved accuracy , consistent with the milder sign problem in the spin-basis formulation. Conclusion —In this w ork, we presen ted a new annealing-based NF metho d that enables reliable sam- pling of the Hubbard mo del in the spin basis at finite c hemical p otential. By learning the sign-quenched the- ory and exploiting the milder sign problem in this for- m ulation, our approac h outp erforms state-of-the-art op- timized HMC calculations b y up to an order of magnitude in b oth accuracy and uncertain ty for system sizes up to N x = 18. The prop osed annealing sc heme enables controlled sampling of multimodal target distributions without in- tro ducing significant training ov erhead or strong dep en- dence on the system size. As a result, the ov erall scala- bilit y of the method is primarily determined b y the prop- erties of the underlying generativ e mo del. While the Re- alNVP architecture employ ed here p erforms well at the lattices sizes considered, its training cost is exp ected to increase for larger volumes. Building on recent work on impro ving the numerical stabilit y of action-based Hub- bard mo del s im ulations [58], we plan to extend our an- nealing sc heme to more expressive and scalable genera- tiv e mo dels in future work. T ak en together, these results establish a broadly ap- plicable framework for efficien t and accurate sampling of the Hubbard mo del in regimes where conv en tional Mon te Carlo metho ds are limited b y ergo dicit y loss or sev ere sign problems. A cknow le dgements —The authors thank Johann Ost- mey er for insightful discussions. The authors gratefully ac knowledge the granted access to the Marvin cluster hosted by the Universit y of Bonn, as well as the comput- ing time granted by the JARA V ergab egremium and pro- vided on the JARA Partition part of the sup ercomputer JURECA at F orsch ungszen trum J ¨ ulic h [59]. This pro ject w as supp orted by the Deutsche F orsch ungsgemeinsc haft (DF G, German Researc h F oundation) as part of the CR C 1639 NuMeriQS – pro ject no. 511713970. ∗ sc huh@hiskp.uni-bonn.de [1] D. P . Arov as, E. Berg, S. A. Kivelson, and S. Raghu, Ann u. Rev. Condens. Matter Phys. 13 , 239 (2022). [2] J. Hubbard, Pro ceedings of the Roy al Society of London. Series A. Mathematical and Physical Sciences 276 , 238 (1963). [3] J. Hubbard, Pro ceedings of the Roy al Society of London. Series A. Mathematical and Physical Sciences 277 , 237 (1964). [4] J. Hubbard, Pro ceedings of the Roy al Society of London. Series A. Mathematical and Physical Sciences 281 , 401 (1964). [5] N. F. Mott and Z. Zinamon, Rep orts on Progress in Ph ysics 33 , 881 (1970). [6] M. Qin, T. Sch¨ afer, S. Andergassen, P . Corb oz, and E. Gull, Ann. Rev. Condensed Matter Ph ys. 13 , 275 (2022), arXiv:2104.00064 [cond-mat.str-el]. [7] B. Xiao, Y.-Y. He, A. Georges, and S. Zhang, Ph ys. Rev. X 13 , 011007 (2023). [8] H. Xu, C.-M. Chung, M. Qin, U. Schollw¨ oc k, S. R. White, and S. Zhang, Science 384 , eadh7691 (2024). [9] L. Razmadze and T. Luu, P oS LA TTICE2024 , 071 (2024). [10] J. Ostmeyer, E. Berko witz, S. Krieg, T. A. L¨ ahde, T. Luu, and C. Urbach, Ph ys. Rev. B 102 , 245105 (2020). [11] J. Ostmeyer, L. Razmadze, E. Berko witz, T. Luu, and U.-G. Meißner, Phys. Rev. B 109 , 195135 (2024). [12] P . Sinilko v, E. Berko witz, T. Luu, and M. Ro dek amp, P oS LA TTICE2024 , 075 (2025). [13] M. Rodek amp, E. Berk owitz, C. G¨ antgen, S. Krieg, T. Luu, J. Ostmeyer, and G. P ederiv a, Eur. Ph ys. J. B 98 (2025), 10.1140/ep jb/s10051-024-00859-1. [14] P . Buivido vic h, D. Smith, M. Ulyb yshev, and L. v on Smek al, Phys. Rev. B 98 , 235129 (2018). [15] S. Beyl, F. Goth, and F. F. Assaad, Phys. Rev. B 97 , 085144 (2018). [16] D. Smith and L. von Smek al, Phys. Rev. B 89 , 195429 (2014). [17] M. V. Ulybyshev, P . V. Buividovic h, M. I. Katsnel- son, and M. I. Polik arp o v, Phys. Rev. Lett. 111 , 056801 (2013). [18] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B 41 , 9301 (1990). [19] M. Rodek amp, E. Berk owitz, C. G¨ antgen, S. Krieg, T. Luu, and J. Ostmeyer, Phys. Rev. B 106 , 125139 (2022), arXiv:2203.00390 [physics.comp-ph]. 6 [20] C. G¨ antgen, E. Berk owitz, T. Luu, J. Ostmeyer, and M. Ro dek amp, Phys. Rev. B 109 , 195158 (2024), arXiv:2307.06785 [cond-mat.str-el]. [21] D. Sc huh, J. Kreit, E. Berko witz, L. F unck e, T. Luu, K. A. Nicoli, and M. Ro dek amp, PoS LA TTICE2024 , 069 (2025). [22] D. Sch uh, J. Kreit, E. Berko witz, L. F unck e, T. Luu, K. A. Nicoli, and M. Rodek amp, (2025), arXiv:2506.17015 [cond-mat.str-el]. [23] J.-L. Wynen, E. Berko witz, C. K¨ orb er, T. A. L¨ ahde, and T. Luu, Ph ys. Rev. B 100 , 075141 (2019), arXiv:1812.09268 [cond-mat.str-el]. [24] R. C. Bro wer, D. Schaic h, and C. Rebbi, PoS Lattice 2011 , 056 (2012). [25] T. Luu and T. A. L¨ ahde, Phys. Rev. B 93 , 155106 (2016). [26] M. Ro dek amp, Machine L e arning for Path Deformation and Bayesian Data A nalysis in Sele cte d Lattic e Field The ories , Ph.d. thesis, Rheinisc he F riedric h-Wilhelms- Univ ersit¨ at Bonn, Bonn (2024). [27] H. F. T rotter, Proceedings of the American Mathematical So ciet y 10 , 545 (1959). [28] M. Suzuki, Progress of Theoretical Physics 56 , 1454 (1976). [29] J. Hubbard, Phys. Rev. Lett. 3 , 77 (1959). [30] P . Buivido vic h, D. Smith, M. Ulyb yshev, and L. v on Smek al, Phys. Rev. B 98 , 235129 (2018). [31] Note that w e could hav e equiv alently defined the chemi- cal potential term on the diagonal of the hopping matrix. [32] S. A. Gottlieb, W. Liu, D. T oussain t, R. L. Renken, and R. L. Sugar, Phys. Rev. D 35 , 2531 (1987). [33] T. Luu and T. A. L¨ ahde, Phys. Rev. B 93 , 155106 (2016). [34] R. T. Scalettar, E. Y. Loh, J. E. Gub ernatis, A. Moreo, S. R. White, D. J. Scalapino, R. L. Sugar, and E. Dagotto, (1989). [35] V. I. Iglovik ov, E. Khatami, and R. T. Scalettar, Phys- ical Review B 92 (2015), 10.1103/ph ysrevb.92.045110. [36] F. T emmen, E. Berko witz, A. Kennedy , T. Luu, J. Ost- mey er, and X. Y u, P oS LA TTICE2024 , 068 (2024). [37] The relation b etw een fermion matrices at finite chemical p oten tial and the emergence of a sign problem is derived in App. D. [38] Z. F o dor and S. D. Katz, Ph ys. Lett. B 534 , 87 (2002), arXiv:hep-lat/0104001. [39] P . de F orcrand, S. Kim, and T. T ak aishi, Nucl. Ph ys. B Pro c. Suppl. 119 , 541 (2003), arXiv:hep-lat/0209126. [40] M. Giordano, K. Kapas, S. D. Katz, D. Nogradi, and A. Pasztor, JHEP 05 , 088 (2020), arXiv:2004.10800 [hep- lat]. [41] Z. F o dor and S. D. Katz, JHEP 03 , 014 (2002), arXiv:hep-lat/0106002. [42] M. T roy er and U.-J. Wiese, Phys. Rev. Lett. 94 , 170201 (2005), arXiv:cond-mat/0408370. [43] G. Aarts, J. Phys. Conf. Ser. 706 , 022004 (2016), arXiv:1512.05145 [hep-lat]. [44] C. Berger, L. Rammelm ¨ uller, A. Loheac, F. Ehmann, J. Braun, and J. Drut, Physics Rep orts 892 , 1 (2021), complex Langevin and other approaches to the sign prob- lem in quantum many-bo dy ph ysics. [45] I. Kobyzev, S. J. Prince, and M. A. Brubak er, IEEE T ransactions on P attern Analysis and Machine Intelli- gence 43 , 3964 (2021). [46] G. Papamak arios, E. Nalisnick, D. J. Rezende, S. Mo- hamed, and B. Lakshminara yanan, J. Mach. Learn. Res. 22 (2021). [47] L. Dinh, D. Krueger, and Y. Bengio, in 3r d Interna- tional Conferenc e on L e arning R epr esentations, ICLR 2015, Workshop T r ack Pr o c e e dings (2015). [48] L. Dinh, J. Sohl-Dic kstein, and S. Bengio, in Interna- tional Confer enc e on L e arning R epr esentations (2017). [49] C. Durk an, A. Bek aso v, I. Murray , and G. P apamak arios, in A dvanc es in Neur al Information Pr o c essing Systems , V ol. 32 (Curran Associates, Inc., 2019). [50] D. P . Kingma and P . Dhariwal, in A dvanc es in Neur al Information Pr o c essing Systems , V ol. 31 (Curran Asso- ciates, Inc., 2018). [51] R. T. Q. Chen, Y. Rubanov a, J. Bettencourt, and D. K. Duv enaud, in A dvanc es in Neur al Information Pr o c essing Systems , V ol. 31 (Curran Asso ciates, Inc., 2018). [52] S. Kullback and R. A. Leibler, The annals of mathemat- ical statistics 22 , 79 (1951). [53] T. Mink a et al. , Diver genc e Me asur es and Message Pass- ing , T echnical Rep ort (Microsoft Research, 2005). [54] D. C. Hack ett, C.-C. Hsieh, S. Pon tula, M. S. Alb ergo, D. Bo yda, J.-W. Chen, K.-F. Chen, K. Cranmer, G. Kan- w ar, and P . E. Shanahan, (2021), [hep-lat]. [55] K. A. Nicoli, C. J. Anders, T. Hartung, K. Jansen, P . Kessel, and S. Nak a jima, Phys. Rev. D 108 , 114501 (2023), arXiv:2302.14082 [hep-lat]. [56] J. Kreit, D. Sch uh, K. A. Nicoli, and L. F unc k e, (2025), arXiv:2505.19619 [cs.LG]. [57] R. M. Neal, Statistics and Computing 11 , 125 (2001). [58] J. Kreit, A. Bulgarelli, L. F unck e, T. Luu, D. Sch uh, S. Singh, and L. V erzic helli, “T o w ard scalable nor- malizing flows for the h ubbard model,” (2026), arXiv:2601.18273 [cond-mat.str-el]. [59] J¨ ulich Sup ercomputing Cen tre, Journal of large-scale re- searc h facilities 7 , A 182 (2021). END MA TTER App endix A: Two-site mo del distribution —The formal ergo dicit y problem in the charge basis formulation of the Hubbard mo del is lifted aw ay from half-filling [23]. Ho w- ev er, this do es not apply to the spin basis, where the ergo cit y problem is of practical nature, i.e., when mo des are widely separated. W e illustrate this on a a tw o-site lattice, where the distribution can b e easily visualized. As shown in Fig. 5, at v anishing chemical p oten tial (left) the distribution takes the form of a correlated Gaus- sian. Introducing a chemical p otential (right) splits the distribution into t wo fairly disconnected mo des, indicat- ing that the practical ergo dicity issues still p ersist. This demonstrates that the main c hallenge faced b y HMC [23] in this form ulation is not remov ed aw a y from half-filling. T o demonstrate that our annealing-based NF metho d can nevertheless mo del this challenging distribution, we pro vide one-b ody correlation functions in Fig. 6, finding excellen t agreement with the exact solution. App endix B: Hybrid Monte Carlo simulations in the spin b asis —HMC is kno wn to struggle with ergo dicity when simulating the Hubbard mo del in the spin basis 7 FIG. 5. Distribution of the tw o-site Hubbard model in the spin basis ( U = 3, β = 5, and N t =20) for µ = 0 (left) and µ = 2 (righ t). Sho wn are the fields summed ov er the temporal direction, i.e., ϕ x = P N t − 1 t =0 ϕ xt . 0.0 0.5 1.0 1.5 t 1 0 2 1 0 1 C ( t ) ED NF SP FIG. 6. One-b o dy correlation functions for the tw o-site Hub- bard mo del with U = 3, β = 5, µ = 2, and N t = 20. Shown are exact diagonalization results at fixed discretization (ED, grey) and those obtained by our annealing-based NF ansatz in the spin basis (NF SP , yello w). [23]. In Fig. 7, we illustrate this for the eight-site setup used in the main text ( U = 2, β = 8, µ = 1 . 75, N t = 40). W e compare HMC to our annealing-based NF ansatz and a ground truth obtained by ED. HMC exhibits the usual non-ergo dic behavior: ov er- sho oting at early times and undersho oting at late times (or vice versa), with v ariations dep ending on initializa- tion. In contrast, the NF metho d pro duces results inde- p enden t of initialization, agreeing with the exact solution and demonstrating ergodic sampling of the target distri- bution. App endix C: F urther r esults —In the main text, we fo- cused on U = 2 for comparison with related work. T o demonstrate applicability at larger U , we sho w results for an eigh t-sites hexagonal lattice with U = 3, β = 6, µ = 2 . 75, and N t = 40 in Fig. 8. Excellen t agreement with the exact solution is observed. App endix D: Non-p ositivity of pr ob ability weight —W e sho w that ev en in the spin basis, a charge-conjugation- lik e symmetry exists at half-filling and is explicitly brok en 0 2 4 6 8 t 3 × 1 0 1 4 × 1 0 1 5 × 1 0 1 6 × 1 0 1 C ( t ) ED HMC SP Run 1 HMC SP Run 2 NF SP Run 1 NF SP Run 2 FIG. 7. One-bo dy correlation functions for the hexagonal eigh t-site Hubbard mo del with U = 2, β = 8, µ = 1 . 75, and N t = 40. Shown are contin uum exact diagonalization results (ED, grey), HMC in the spin basis with tw o initial- izations (HMC SP Run 1/2, blue/grey), and our annealing- based NF ansatz with tw o initializations (NF SP Run 1/2, y ellow/orange). Small deviations from the exact results at early times are due to discretization errors. 0 1 2 3 4 5 6 t 1 0 3 1 0 2 1 0 1 1 0 0 C ( t ) ED NF SP FIG. 8. One-bo dy correlation functions for the hexagonal eigh t-site Hubbard mo del with U = 3, β = 6, µ = 2 . 75, and N t = 40. Shown are contin uum exact diagonalization results (ED, grey) and 1M samples generated b y our anneal ing-based NF ansatz (NF SP , yello w). at finite c hemical potential, revealing the origin of the sign problem in this formulation. Extending the treatment of the fermion determinant presen ted in Ref. [23] to finite c hemical potential, w e first consider the case ˜ µ = 0: det M [ ϕ | ˜ κ ] ≡ det M [ ϕ | ˜ κ, ˜ µ = 0] . (15) In the exp onential discretization, the fermion determi- nan t can b e written as det M [ ϕ | ˜ κ ] = det( I N x + A ) , (16) where A ≡ T N t − 1 T N t − 2 · · · T 0 and T t ≡ e h F t with F t [ ϕ ] x ′ x = e ϕ x,t − 1 δ x ′ x . Equiv alen tly , it can b e expressed 8 as det M [ ϕ | ˜ κ ] = det( A ) · det  A − 1 + I N x  . (17) Here, det e h = 1 since h is traceless, and the pro duct ov er the field exp onentials giv es det( F N t − 1 · · · F 0 ) = Y t ′ e P x ϕ x,t ′ = e P x,t ϕ x,t = e Φ , (18) leading to det M [ ϕ | ˜ κ ] = e Φ det  I N x + A − 1  . (19) Using F − 1 t [ ϕ ] = F t [ − ϕ ], the bipartite condition det M [ ϕ | ˜ κ ] = det M [ ϕ | − ˜ κ ] [23], and the in v ariance det( I + X ) = det  I + X T  for any square matrix X , we can rearrange the factors in A . Since F t [ ϕ ] is diagonal and h is symmetric, we hav e F t [ ϕ ] T = F t [ ϕ ] and ( e h ) T = e h , giving det M [ − ϕ | ˜ κ ] = det  I N x + A − 1  . (20) Comparing Eqs. (19) and (20) yields det M [ ϕ | ˜ κ ] = e Φ det M [ − ϕ | ˜ κ ] , (21) ensuring p ositivit y of the probability weigh t in the spin basis: W [ ϕ | ˜ κ ] = det M [ ϕ | ˜ κ ] det M [ − ϕ | ˜ κ ] e − 1 2 ˜ U P x,t ϕ 2 x,t = e Φ det( M [ − ϕ, ˜ κ ]) 2 e − 1 2 ˜ U P x,t ϕ 2 x,t > 0 . (22) A t finite chemical p otential, following Ref. [20], w e mo d- ify F t → ˜ F t [ ϕ, ˜ µ ] x ′ x = e ϕ x, ( t − 1) + ˜ µ δ x ′ x . (23) With this, and using Eqs. (18) and (19), we obtain det M [ ϕ | ˜ κ, ˜ µ ] = e Φ+ V ˜ µ det  I N x + ˜ A − 1  , (24) where V = N t × N x . Aside from the harmless factor e V ˜ µ , the only mo dification app ears in ˜ A − 1 , constructed from ˜ F t [ − ϕ, − ˜ µ ]. The loss of p ositivity originates from ˜ A − 1 , introducing factors of e − ˜ µ . Applying the bipartite condition yields det M [ ϕ | ˜ κ, ˜ µ ] = e Φ+ V ˜ µ det M [ − ϕ | ˜ κ, − ˜ µ ] . (25) The main obstacle to p ositivity is that Eq. (25) relates a determinan t at p ositive c hemical p oten tial to one at neg- ativ e c hemical p oten tial. Since b oth spin sp ecies couple to the same chemical p otential, the resulting probability w eight W [ ϕ | ˜ κ, ˜ µ ] = det M [ ϕ | ˜ κ, ˜ µ ] det M [ − ϕ | ˜ κ, ˜ µ ] e − 1 2 ˜ U P x,t ϕ 2 x,t = e Φ+ V ˜ µ det  M [ − ϕ | ˜ κ, − ˜ µ ] M [ − ϕ | ˜ κ, ˜ µ ]  × e − 1 2 ˜ U P x,t ϕ 2 x,t , (26) no longer factorizes into a manifestly p ositive form, ex- plicitly demonstrating the origin of the sign problem.