Integral formula for the propagator of the one-dimensional Hubbard model

In tegral form ula for the propagator of the one-dimensional Hubbard mo del T aiki Ishiy ama 1 , Kazuy a F ujimoto 1 , and T omohiro Sasamoto 1 1 Dep artment of Physics, Institute of Scienc e T okyo, 2-12-1 Ookayama, Me gur o-ku, T okyo 152-8551, Jap an F ebruary 26, 2026 Abstract W e present an exact integral form ula for the m ulti-particle propagator of the one- dimensional F ermi–Hubbard mo del on an infinite lattice. The pro of is based on the nested Bethe ansatz without relying on the string hypothesis. Our form ula enables an explicit in tegral representation of the time evolution of arbitrary finite- particle wa ve functions and thereby pro vides a foundation for the exact analysis of nonequilibrium dynamics in the Hubbard mo del. It can further b e applied to related open quantum mo dels. 1 In tro duction The one-dimensional (1D) F ermi–Hubbard mo del serves as a fundamental platform for studying strongly correlated electron systems in one dimension, where quan tum fluctua- tions enhanced b y lo w dimensionality giv e rise to in triguing phenomena such as T omonaga– Luttinger liquid b eha vior and spin-charge separation [1]. A distinctive feature of this mo del is its in tegrabilit y . Lev eraging this prop erty , since the seminal work of Lieb and W u [2], its spectral prop erties and thermo dynamics hav e b een extensiv ely in vestigated using the nested Bethe ansatz. In recen t y ears, interest in the 1D F ermi–Hubbard mo del has shifted significan tly from stationary prop erties to nonequilibrium dynamics, motiv ated b y remark able exp erimen- tal and theoretical adv ancements. On the exp erimental side, state-of-the-art ultracold atom experiments hav e enabled the faithful sim ulation of real-time ev olution in isolated quan tum many-bo dy systems [3, 4]. In particular, Hubbard-type mo dels can b e naturally realized in such platforms, where the dimensionality , particle statistics, and in teraction strength can b e widely tuned. On the theoretical side, tensor-netw ork simulations (TNS), including the time-evolving blo c k decimation metho d [5] and time-dep endent densit y ma- trix renormalization group [6], ha ve emerged as efficien t n umerical to ols for sim ulating the time ev olution of 1D quantum many-bo dy systems [7, 8]. F urthermore, generalized h y- dro dynamics (GHD) [9 – 11] has b een formulated, whic h describ es hydrodynamic b ehavior at Euler scales for general integrable systems. While these theoretical framew orks ha ve successfully elucidated nonequilibrium prop- erties of the 1D F ermi–Hubbard mo del [12 – 22], w e note that the accessible time scale of 1 TNS is strongly limited by entanglemen t growth, and GHD primarily describ es the av er- age ballistic transp ort of conserved c harges. Consequently , exact microscopic approaches are indisp ensable to complemen t these framew orks, particularly to capture phenomena b ey ond the regimes of v alidity of these metho ds. Despite the in tegrabilit y of the 1D F ermi–Hubbard mo del, p erforming exact micro- scopic calculations of nonequilibrium dynamics remains a formidable challenge. F or sim- pler quantum integrable systems, suc h as the Lieb–Liniger Bose gas and the XXZ spin c hain, exact in tegral representations of the m ulti-particle propagator on an infinite in- terv al, often referred to as the Y udson represen tation [23, 24], are well established via the Bethe ansatz [25 – 29]. These form ulas enable an exact analysis of time evolution for finite particles on the infinite interv al without relying on the string h yp othesis. More recen tly , several studies hav e attempted to in vestigate the nonequilibrium dynamics of finite-densit y initial states in the XXZ spin c hain b y taking the thermodynamic limit of these integral represen tations [30, 31]. In parallel, closely related integral form ulas for m ulti-particle propagators hav e b een developed indep enden tly in the field of integrable sto c hastic interacting systems [32 – 35]. In this context, v arious studies hav e demonstrated that the use of such form ulas together with the Marko v dualit y [36] enables the anal- ysis of nonequilibrium dynamics at finite density [37 – 41]. Ho w ever, no exact integral represen tation of the propagator is known for the 1D F ermi–Hubbard mo del. In this work, w e derive an exact integral form ula for the multi-particle propagator of the 1D F ermi–Hubbard mo del on an infinite lattice. The formula is expressed as a m ultiple con tour in tegral of the Bethe wa ve function with appropriately c hosen in tegra- tion contours. Unlik e the Lieb–Liniger mo del and the XXZ spin c hain, the presence of in ternal degrees of freedom in the 1D F ermi–Hubbard mo del necessitates the use of the nested Bethe ansatz [1, 2]. Consequently , the deriv ation of the formula is technically more in v olved, requiring us to explicitly emplo y the nested structure of the Bethe w av e func- tion. The form ula pro vides an explicit representation of the time ev olution of arbitrary finite-particle w a v e functions on the infinite lattice. It can furthermore b e applied to op en quan tum systems whose effectiv e description is g iven by the Hubbard mo del with complex in teraction [42 – 45]. Th us our formula establishes a foundation for the exact microscopic analysis of nonequilibrium dynamics in these systems. The pap er is organized as follows. In section 2, we define the mo del and the multi- particle propagator. In section 3, we briefly review the nested Bethe ansatz for the 1D F ermi–Hubbard mo del [1, 2]. In section 4, w e deriv e an in tegral form ula for the propagator using the nested Bethe ansatz [1, 2]. W e conclude the pap er in section 5. 2 Setup: the 1D F ermi–Hubbard mo del In this work, we consider the 1D F ermi–Hubbard mo del on an infinite lattice Z . F or brevit y , we refer to it simply as the Hubbard mo del throughout the rest of this paper. The Hamiltonian is giv en by ˆ H := − X x ∈ Z X σ = ↓ , ↑  ˆ c † x,σ ˆ c x +1 ,σ + ˆ c † x +1 ,σ ˆ c x,σ  + 4 u X x ∈ Z ˆ n x, ↓ ˆ n x, ↑ . (1) Here, ˆ c x,σ (ˆ c † x,σ ) and ˆ n x,σ := ˆ c † x,σ ˆ c x,σ denote the annihilation (creation) op erator and the n um b er op erator of fermions at site x with spin σ , resp ectively . The on-site in teraction strength is parameterized b y 4 u , where the factor 4 is introduced for later conv enience. While the in teraction parameter u is usually real in (1), w e allo w it to tak e complex v alues 2 u ∈ C \ { 0 } throughout this work. W e note that the non-in teracting case is excluded due to technical reasons that will b ecome clear in section 4. The Hamiltonian conserves the total particle num b er ˆ N and the num b er of down-spin particles ˆ N ↓ , [ ˆ H , ˆ N ] = [ ˆ H , ˆ N ↓ ] = [ ˆ N , ˆ N ↓ ] = 0 , (2) where ˆ N := P x ∈ Z P σ = ↓ , ↑ ˆ n x,σ and ˆ N ↓ := P x ∈ Z ˆ n x, ↓ . Hereafter, w e fo cus on the sector with N particles and M down-spin particles, assuming N , M > 0. The quantit y of interest is the propagator ψ t ( x ; a | y ; b ) := ⟨ x ; a | e − it ˆ H | y ; b ⟩ , (3) where we define the W annier state with the v acuum state | 0 ⟩ as | x ; a ⟩ := ˆ c † x N ,a N · · · ˆ c † x 1 ,a 1 | 0 ⟩ . (4) The propagator describ es the transition amplitude that the system, initially prepared in the state | y ; b ⟩ , ev olves into the state | x ; a ⟩ at time t . In terms of the propagator, any N -particle state | Φ( t ) ⟩ can b e expanded as follows, | Φ( t ) ⟩ = 1 ( N !) 2 X x , y ∈ Z N X a , b ∈{↓ , ↑} N | x ; a ⟩⟨ x ; a | e − it ˆ H | y ; b ⟩⟨ y ; b | Φ(0) ⟩ (5) = 1 ( N !) 2 X x , y ∈ Z N X a , b ∈{↓ , ↑} N ψ t ( x ; a | y ; b ) ⟨ y ; b | Φ(0) ⟩| x ; a ⟩ . (6) Namely , one can exactly calculate any N -particle state | Φ( t ) ⟩ if one obtains the exact solution of ψ t ( x ; a | y ; b ). The propagator satisfies the following initial condition, ψ t ( x ; a | y ; b ) | t =0 = det[ δ x j ,y k δ a j ,b k ] N j,k =1 , (7) and the Sc hr¨ odinger equation of the one-dimensional Hubbard mo del in first quantiza- tion [1], i d dt ψ t ( x ; a | y ; b ) = H N ψ t ( x ; a | y ; b ) , (8) H N := − N X j =1 (∆ + j + ∆ − j ) + 4 u X 1 ≤ j 1, assuming that (30) holds for M − 1. F or n = 1 , . . . , M , define the subset of S M as S ( n ) M := { R ∈ S M | R (1) = n } . (36) The sets S ( n ) M ( n = 1 , . . . , M ) are disjoin t and their union giv es the symmetric group, S M = M G n =1 S ( n ) M . (37) This decomp osition allo ws us to ev aluate the total sum ov er S M b y separately considering con tributions from each S ( n ) M . No w consider I s ( R ) for R ∈ S (1) M . In this case, the λ 1 -in tegral can b e ev aluated as in the case M = 1. After the λ 1 -in tegration, we hav e I s ( R ) = δ α 1 ,β 1 M Y n =2 " I Γ ( n ) s dλ n 2 π i 2 iu ( λ n − s β n − iu )( λ n − s α R − 1 ( n ) + iu ) # × Y 2 ≤ m 1). In this case, we will show the following iden tit y , X R ∈ S ( n ) M I s ( R ) = 0 for n > 1 . (39) 10 When β n ≤ α 1 , it follows that β 1 < β n ≤ α 1 < α R − 1 (1) . (40) Therefore, the integrand in (31) is holomorphic inside the λ 1 -con tour, which prov es (39). Henceforth, we assume α 1 < β n and prov e (39) in this case. F or R ∈ S ( n ) M , I s ( R ) can b e expressed as I s ( R ) = I Γ ( n ) s dλ n 2 π i 2 iu ( λ n − s β n − iu )( λ n − s α 1 + iu ) n − 1 Y m =1 λ m − λ n + 2 iu λ m − λ n − 2 iu β n − 1 Y β = α 1 λ n − s β + iu λ n − s β − iu × ( λ n -indep enden t part) . Based on (35), w e ev aluate the λ n -in tegral b y pic king up the residues at λ n = λ m − 2 iu ( m = 1 , . . . , n − 1), which leads to I s ( R ) = n − 1 X m =1 M Y k =1 k  = n h I Γ ( k ) s dλ k 2 π i 2 iu ( λ k − s β k − iu )( λ k − s α R − 1 ( k ) + iu ) i × − 8 u 2 ( λ m − s β n − 3 iu )( λ m − s α 1 − iu ) n − 1 Y l =1 l  = m λ l − λ m + 4 iu λ l − λ m β n − 1 Y β = α 1 λ m − s β − iu λ m − s β − 3 iu × Y 1 R − 1 ( l ) λ k − λ l + 2 iu λ k − λ l − 2 iu × M Y k =1 k  = n h α R − 1 ( k ) − 1 Y α =1 λ k − s α − iu λ k − s α + iu β k − 1 Y β =1 λ k − s β + iu λ k − s β − iu i . W e next p erform the integration with resp ect to λ m . F rom the conditions α 1 < α R − 1 ( m ) and β m < β n , the integrand of the ab ov e expression has no singularities at λ m = s j + iu ( j = 1 , . . . , N ). F urthermore, giv en our c hoice of the contour and the condition s j − s k  = 2 iu for all j, k , no p oles other than λ m = λ l ( l = 1 , . . . , n − 1; l  = m ) can b e lo cated inside Γ ( m ) s . Hence, the only p ossible singularities within the λ m -con tour are the p oints λ m = λ l . Thus, after the λ m -in tegration, one obtains I s ( R ) = n − 1 X m =1 n − 1 X l =1 l  = m T s ( l , m ; R ) (41) where we define T s ( l , m ; R ) := M Y k =1 k  = n,m h I Γ ( k ) s dλ k 2 π i 2 iu ( λ k − s β k − iu )( λ k − s α R − 1 ( k ) + iu ) i × 2 iu χ [ λ l ∈ In t Γ ( m ) s ] ( λ l − s β m − iu )( λ l − s α R − 1 ( m ) + iu ) × Res λ m = λ l " − 8 u 2 ( λ m − s β n − 3 iu )( λ m − s α 1 − iu ) n − 1 Y k =1 k  = m λ k − λ m + 4 iu λ k − λ m β n − 1 Y β = α 1 λ m − s β − iu λ m − s β − 3 iu × Y 1 R − 1 ( k ) λ j − λ k + 2 iu λ j − λ k − 2 iu × M Y k =1 k  = n h α R − 1 ( k ) − 1 Y α =1 λ k − s α − iu λ k − s α + iu β k − 1 Y β =1 λ k − s β + iu λ k − s β − iu i # (42) 11 with the indicator function χ [ · ] and the interior In t Γ ( m ) s of the λ m -con tour. No w we observe that T s ( l , m ; R ) + T s ( l , m ; Π l,m R ) = 0 . (43) This relation follows from the fact that, up on taking the residue at λ m = λ l , the first factor in the fourth line of (42) for Π l,m R is the negativ e of that for R , whereas all the other factors coincide. Therefore, using the ab o ve relation and the decomp osition of S ( n ) M b y even and o dd p ermutations as S ( n ) M = S ( n ) M , even F Π l,m S ( n ) M , even , we obtain the desired iden tit y (39) as follo ws, X R ∈ S ( n ) M I s ( R ) = X R ∈ S ( n ) M n − 1 X m =1 n − 1 X l =1 l  = m T s ( l , m ; R ) = n − 1 X m =1 n − 1 X l =1 l  = m X R ∈ S ( n ) M , even [ T s ( l , m ; R ) + T s ( l , m ; Π l,m R )] = 0 . Finally , combining (38) and (39), we obtain the identit y (30), whic h completes the pro of of Prop osition 1 and thereby establishes Lemma 1. 4.2 Pro of of Lemma 2 W e complete the pro of of Step (b) by proving Lemma 2. F or notational conv enience, w e denote the summand on the left-hand side of (25) as I c ( P ). Then the iden tit y (25) can b e written as X P ∈ S N I c ( P ) = det[ δ x j ,y k δ a j ,b k ] N j,k =1 . (44) As sho wn b elow, Lemma 2 can b e pro ved in a manner similar to the proof of Propo- sition 1. This is b ecause the amplitude A ( s ) P ( a | b ) giv en in (24) is comp osed of the Y - op erators (18), Y ( s j − s k ) = ( s j − s k )Π + 2 iu s j − s k + 2 iu , (45) whose structure is analogous to that of the factors ( λ j − λ k + 2 iu ) / ( λ j − λ k − 2 iu ) appearing in (31). In the presen t case, ho wev er, these factors are matrices acting on the auxiliary spin space, which makes the calculation more in volv ed. T o handle this complication, it is useful to in tro duce the graphical represen tation of A ( s ) P ( a | b ). The expression (24) admits a natural graphical represen tation, whic h pro vides an intuitiv e understanding of the structure of the amplitude A ( s ) P ( a | b ). Figure 2 shows the corresp onding diagrams for (a) P = Π 1 , 2 and (b) P = Π 1 , 3 . Eac h arrow is lab eled with the rapidit y s j , and the p ositions where it starts and ends corresp ond to site j and site P − 1 ( j ). Eac h in tersection indicates the position where the Y -op erator acts and whic h pair of rapidities is scattered by the Y -op erator. F rom the sequence and arrangement of these in tersections, one can directly reconstruct the corresp onding amplitude. Note that 12 Figure 2: Graphical represen tations of A ( s ) P ( a | b ) for (a) P = Π 1 , 2 . (b) P = Π 1 , 3 . F rom Lemma 1, one has A ( s ) Π 1 , 2 ( a | b ) = ⟨ a | Y 1 , 2 ( s 1 − s 2 ) | b ⟩ and A ( s ) Π 1 , 3 ( a | b ) = ⟨ a | Y 1 , 2 ( s 2 − s 3 ) Y 2 , 3 ( s 1 − s 3 ) Y 1 , 2 ( s 1 − s 2 ) | b ⟩ . the (braid) Y ang–Baxter relation (26) guaran tees the in v ariance of the amplitude under a change in the order of intersections. W e pro ceed by induction on N . The case N = 1 is obvious. Let us consider N = 2. In this case, one has I c (id) = δ x , y δ a , b straigh tforw ardly . F or I c (Π 1 , 2 ), we make the substitution z 2 → z 2 /z 1 , which leads to I c (Π 1 , 2 ) = − I | z 1 | = r dz 1 2 π z 1 I | z 2 | = r dz 2 2 π iz 2 ⟨ a | ( z 2 1 − 1 − z 2 + z 2 1 /z 2 )Π 1 , 2 − 4 uz 1 z 2 1 − 1 − z 2 + z 2 1 /z 2 − 4 uz 1 | b ⟩ z x 1 − y 2 2 z x 2 − x 1 + y 2 − y 1 1 . Since | z 2 | = r on the z 2 -con tour, 1 / ( z 2 1 − 1 − z 2 + z 2 1 /z 2 − 4 uz 1 ) is analytic inside the z 1 - con tour. The term z x 2 − x 1 + y 2 − y 1 1 is also analytic from the conditions x 1 ≤ x 2 and y 1 ≤ y 2 . Therefore, from the residue theorem, it follo ws that I c (Π 1 , 2 ) = − I | z 2 | = r dz 2 2 π iz 2 z x 1 − y 2 2 δ a 1 ,b 2 δ a 2 ,b 1 δ x 1 ,x 2 δ y 1 ,y 2 (46) = − δ x 1 ,y 2 δ x 2 ,y 1 δ y 1 ,y 2 δ a 1 ,b 2 δ a 2 ,b 1 (47) = − δ x 1 ,y 2 δ x 2 ,y 1 δ a 1 ,b 2 δ a 2 ,b 1 . (48) In the last line, w e use the iden tity δ x 1 ,y 2 δ x 2 ,y 1 δ y 1 ,y 2 = δ x 1 ,y 2 δ x 2 ,y 1 whic h holds when x 1 ≤ x 2 and y 1 ≤ y 2 . Thus one obtains the iden tity (44) for the case N = 2. W e shall pro ve (44) for N ≥ 3 assuming that it holds for N − 1 and N − 2. As in the case of Prop osition 1, we use the decomp osition of the symmetric group: S N = F N n =1 S ( n ) N . The definition of S ( n ) N is given in (36). In what follo ws, w e establish the follo wing iden tities, X P ∈ S (1) N ∪ S (2) N I c ( P ) = det[ δ x j ,y k δ a j ,b k ] N j,k =1 , (49) X P ∈ S ( n ) N I c ( P ) = 0 for n ≥ 3 . (50) These tw o identities prov e (44) for N ≥ 3. The following expression for the amplitude with P ∈ S ( n ) N pla ys an imp ortan t role in pro ving (49) and (50): A ( s ) P ( a | b ) = ⟨ a | · · · − → Y 1 ≤ m 2, it follo ws that x P − 1 (1) − x 1 + y 2 − y 1 − 1 ≥ 0. Therefore, from the residue theorem, we hav e I c ( P ) = 0 when P − 1 (1) > 2 . (53) On the other hand, for the case P − 1 (1) = 2, I c ( P ) can b e ev aluated as I c ( P ) = − sgn( ˜ P ) δ x 1 ,y 2 δ x 2 ,y 1 δ a 1 ,b 2 δ a 2 ,b 1 × h N − 2 Y j =1 I | z j | = r N − j dz j 2 π iz j i A ( s ) ˜ P ( a 3 , · · · , a N | b 3 , · · · , b N ) N − 2 Y j =1 z x j +2 − y ˜ P ( j )+2 ˜ P ( j ) , where w e define ˜ P ∈ S N − 2 as ˜ P ( j ) := P ( j + 2) − 2, ( j = 1 , . . . , N − 2), and used the iden tit y δ x 1 ,y 2 δ x 1 ,x 2 δ y 1 ,y 2 = δ x 1 ,y 2 δ x 2 ,y 1 for x 1 ≤ x 2 and y 1 ≤ y 2 . Therefore, from the induction hypothesis, it follo ws that X P ∈ S (2) N I ( P ) = − δ x 1 ,y 2 δ x 2 ,y 1 δ a 1 ,b 2 δ a 2 ,b 1 det[ δ x j ,y k δ a j ,b k ] N j,k =3 . (54) 14 Com bining (52) and (54), we ha v e X P ∈ S (1) N ,S (2) N I c ( P ) = δ x 1 ,y 1 δ a 1 ,b 1 det[ δ x j ,y k δ a j ,b k ] N j,k =2 − δ x 1 ,y 2 δ x 2 ,y 1 δ a 1 ,b 2 δ a 2 ,b 1 det[ δ x j ,y k δ a j ,b k ] N j,k =3 . On the other hand, the determinant can b e expanded as follows, det[ δ x j ,y k δ a j ,b k ] N j,k =1 = X P ∈ S N sgn( P ) δ x 1 ,y P (1) · · · δ x N ,y P ( N ) δ a 1 ,b P (1) · · · δ a N ,b P ( N ) = δ x 1 ,y 1 δ a 1 ,b 1 det[ δ x j ,y k δ a j ,b k ] N j,k =2 − δ x 1 ,y 2 δ x 2 ,y 1 δ a 1 ,b 2 δ a 2 ,b 1 det[ δ x j ,y k δ a j ,b k ] N j,k =3 . In the last line, we used the fact that Q N j =1 δ x j ,y P ( j ) δ a j ,b P ( j ) v anishes unless P (1) = 1, or P (1) = 2 and P (2) = 1. This restriction arises from the ordering x 1 ≤ · · · ≤ x N , y 1 ≤ · · · ≤ y N , and the distinctness of the pairs ( x j , a j ) and ( y j , b j ). Th us we obtain (49). W e next prov e (50). Let us consider I c ( P ) for P ∈ S ( n ) N , ( n ≥ 3). When x 2 > y 1 , one can verify that I c ( P ) = 0 as the integrand is analytic inside the z 1 -con tour. Hence, in the follo wing, we assume y 1 − x 2 ≥ 0. W e first ev aluate the z n -in tegral in I c ( P ). After the substitution z n → − 1 /z n , I c ( P ) b ecomes I c ( P ) = ( − 1) x 1 − y n sgn( P ) h N Y j =1 j  = n I | z j | = r N − j dz j 2 π iz j i I | z n | =1 /r N − n dz n 2 π i A ( s ) P ( a | b ) z y n − x 1 − 1 n N Y j =1 j  = n z x P − 1 ( j ) − y j j Note that since the function s n = ( z n − 1 /z n ) / 2 i is in v arian t under the substitution, A ( s ) P ( a | b ) also remains unchanged. F rom the assumption y 1 − x 2 ≥ 0, it follows that y n − x 1 − 1 ≥ 0 for n ≥ 3. Hence the z n -singularit y of the in tegrand originates solely from A ( s ) P ( a | b ). In (51), the z n -dep enden t factor takes the following explicit form, Y m,m +1 ( s m − s n ) = − 2 iz n 2 iu + ( s m − s n )Π m,m +1 z 2 n − 2 i ( s m + 2 iu ) z n − 1 . (55) The p oles of this expression are lo cated at z n = ρ ( z m ) and z n = − 1 /ρ ( z m ), where ρ ( z ) is the solution to ρ 2 − 2 i ( s ( z ) + 2 iu ) ρ − 1 = 0 satisfying ρ ( z ) ≃ z for | z | ≪ 1. The p ole z n = ρ ( z m ) lies inside the z n -con tour, while the other lies outside. Hence, by p erforming the integration with resp ect to z n , we obtain I c ( P ) = ( − 1) x 1 − y n sgn( P ) n − 1 X m =1 h N Y j =1 j  = n,m I | z j | = r N − j dz j 2 π iz j i I | z m | = r N − m dz m 2 π i 2 i ρ 2 m − 1 × Res s n = s m +2 iu [ A ( s ) P ( a | b )] z x P − 1 ( m ) − y m − 1 m ρ y n − x 1 +1 m N Y j =1 j  = n,m z x P − 1 ( j ) − y j j , (56) where we defined the shorthand notation ρ m := ρ ( z m ) and used the relation Res z n = ρ m Y m,m +1 ( s m − s n ) = 2 iρ 2 m ρ 2 m − 1 Res s n = s m +2 iu Y m,m +1 ( s m − s n ) . (57) 15 W e next p erform the z m -in tegration in (56). Since ρ m = z m + O ( z 2 m ) for | z m | ≪ 1, the following quan tity , 1 ρ 2 m − 1 z x P − 1 ( m ) − y m − 1 m ρ y n − x 1 +1 m , (58) is holomorphic inside the z m -con tour. Hence again the singularity of the integrand origi- nates solely from Res s n = s m +2 iu A ( s ) P ( a | b ). The follo wing factor in this quantitiy ma y give rise to p oles inside the z m -con tour, ( s l − s m )Π + 2 iu s l − s m + 2 iu (1 ≤ l < m and P − 1 ( m ) < P − 1 ( l )) , (59) ( s l − s n )Π + 2 iu s l − s n + 2 iu    s n = s m +2 iu = ( s l − s m − 2 iu )Π + 2 iu s l − s m , (1 ≤ l < m ) . (60) See figure 3 for the sc hematic illustration of these factors. Note that, with our choice of the contours, the following z m -dep enden t factors are analytic inside the z m -con tour, ( s m − s j )Π + 2 iu s m − s j + 2 iu ( m < j ) . (61) Figure 3: Sc hematic illustration of Res s n = s m +2 iu [ A ( s ) P ( a | b )]. The blue circles and diamonds represen t the factors giv en in (59) and (60), resp ectively , while the red square corresponds to the p ole-con tributing factor in the z n -in tegration. The dashed triangle represen ts the factor given in (62). In terestingly the denominator of (59) cancels up on combination with other factors in Res s n = s m +2 iu A ( s ) P ( a | b ). In other words, we do not need to account for the corresp onding p oles in the in tegration. The k ey p oin t is that, by virtue of the Y ang-Baxter relation (26), one can choose a decomposition of the p ermutation P into adjacent transp ositions suc h that the following sequence app ears in Res s n = s m +2 iu A ( s ) P ( a | b ): Y j,j +1 ( s l − s m ) Y j − 1 ,j ( s l − s m − 2 iu )[ − 2 iu (1 − Π j,j +1 )] (62) = ( s l − s m )Π j,j +1 + 2 iu s l − s m + 2 iu ( s l − s m − 2 iu )Π j − 1 ,j + 2 iu s l − s m [ − 2 iu (1 − Π j,j +1 )] = ( s l − s m + 2 iu − 2 iu )Π j,j +1 + 2 iu s l − s m + 2 iu ( s l − s m + 2 iu − 4 iu )Π j − 1 ,j + 2 iu s l − s m [ − 2 iu (1 − Π j,j +1 )] where j is some integer, and the remaining factors are regular at s m = s l + 2 iu . See app endix C for an explicit example of suc h a decomposition. This sequence corresp onds 16 to a triangle formed by the lines l , m , and n in the graphical representation (see figure 3). F rom the final expression ab ov e, one sees that if (1 − Π j,j +1 )(2Π j − 1 ,j − 1)(1 − Π j,j +1 ) = 0 , (63) then Res s n = s m +2 iu A ( s ) P ( a | b ) is regular at s m = s l + 2 iu . In fact, (63) can b e sho wn by the direct ev aluation of its action on arbitrary basis vectors | a j − 1 , a j , a j +1 ⟩ . Consequently , in the z m -in tegration in (56), the only contributing poles are those of (60). Therefore, we obtain I c ( P ) = n − 1 X m =1 m − 1 X l =1 T c ( l , m ; P ) , (64) where T c ( l , m ; P ) := ( − 1) x 1 − y n sgn( P ) h N Y j =1 j  = n,m I | z j | = r N − j dz j 2 π iz j i 2 i ρ 2 l − 1 2 i z 2 l − 1 × Res s m = s l  Res s n = s m +2 iu [ A ( s ) P ( a | b )]  z x P − 1 ( l ) + x P − 1 ( m ) − y l − y m +1 l ρ y n − x 1 +1 l N Y j =1 j  = n,m,l z x P − 1 ( j ) − y j j . No w we observe that T c ( l , m ; P ) + T c ( l , m ; Π l,m P ) = 0 , (65) whic h is analogous to (43). T o pro v e this, it suffices to show that Res s m = s l  Res s n = s m +2 iu [ A ( s ) P ( a | b )]  = Res s m = s l  Res s n = s m +2 iu [ A ( s ) Π l,m P ( a | b )]  . (66) Although the v alidit y of (66) can b e seen from figure 4, we show it explicitly using the Y -op erator in app endix C. Figure 4: Schematic illustration of (66). The red square and blue diamonds represent the p ole-con tributing factors in the z n - and z m -in tegrations, resp ectively . The intersection b et ween the lines l and m corresp onds to the Y -op erator Y ( s l − s m ). Note that this op erator reduces to the identit y op erator when s l = s m . 17 Using (65) and the decomp osition of S ( n ) N b y even and o dd p ermutations as S ( n ) N = S ( n ) N , even F Π l,m S ( n ) N , even , we obtain X P ∈ S ( n ) N I c ( P ) = X P ∈ S ( n ) N n − 1 X m =1 m − 1 X l =1 T c ( l , m ; P ) = n − 1 X m =1 m − 1 X l =1 X P ∈ S ( n ) N , even [ T c ( l , m ; P ) + T c ( l , m ; Π l,m P )] = 0 Th us, we arriv e at the identit y (50), which, together with the iden tity (49), completes the pro of of Lemma 2. Finally , from Lemma 1 and Lemma 2, w e establish the pro of of Step (b), thereby completing the pro of of Theorem 1. 5 Conclusion In this w ork, w e ha v e presented an exact integral form ula for the multi-particle propa- gator of the one-dimensional F ermi–Hubbard mo del on the infinite lattice. Through this form ula, the time dep endence of arbitrary finite-particle wa ve functions can b e computed explicitly in the Hubbard mo del. The formula can also b e applied to op en quan tum mo dels related to the Hubbard mo del [42 – 45]. Thus, it establishes a microscopic starting p oin t for the exact analysis of nonequilibrium dynamics in these systems. The pro of of the formula has been established by verifying tw o conditions: the satis- faction of the Sc hr¨ odinger equation and the fulfillmen t of the initial condition. While the former follo ws directly from the nested Bethe ansatz [1, 2], the latter presen ts a non-trivial c hallenge. W e resolved this b y lev eraging the nested structure of the Bethe w av e function; sp ecifically , w e p erformed successive residue calculations with resp ect to b oth spin and c harge rapidities. Crucially , our deriv ation works directly on the infinite lattice, thereby a v oiding reliance on the string hypothesis. As an outlook, it would b e in triguing to inv estigate the nonequilibrium dynamics of the Hubbard mo del and related op en quantum systems using the present formalism. Another natural direction is to extend the presen t approach to the Maassarani mo del [49], whic h is the SU( n ) generalization of the Hubbard mo del. Ac knowledgmen ts TI, KF, and TS are grateful to F abian Essler, Cristian Giardin` a, T ak ashi Imam ura, Matteo Macc hini for helpful discussions and commen ts. The w ork of TI has b een supp orted b y JST SPRING, Japan Grant Number JPMJSP2180. The w ork of KF has b een supp orted b y JSPS KAKENHI Gran t No. JP23K13029. The work of TS has b een supp orted by JSPS KAKENHI Grants No. JP21H04432, No. JP22H01143, and No. JP23K22414. 18 A Op en quantum mo dels related to the Hubbard mo del W e explain op en quan tum mo dels relev an t to the Hubbard mo del with complex in ter- action strength and the application of the propagator (3) to these mo dels. In this ap- p endix, w e consider op en quantum systems whose time evolution is gov erned by the Gorini-Kossak o wski-Sudarshan-Lindblad (GKSL) equation [50, 51]: d dt ˆ ρ ( t ) = L [ ˆ ρ ( t )] := − i [ ˆ H , ˆ ρ ( t )] + X x h 2 ˆ L x ˆ ρ ( t ) ˆ L † x − { ˆ L † x ˆ L x , ˆ ρ ( t ) } i . (67) Here, L is the sup erop erator referred to as the Liouvillian, ˆ H denotes the system’s Hamil- tonian, and the Lindblad op erator ˆ L x represen ts the dissipation caused by coupling to the en vironmen t. F or a comprehensive introduction to op en quantum systems, w e refer the reader to [52]. A.1 Tigh t-binding c hain with dephasing noise A tight-binding chain with dephasing noise is describ ed b y the GKSL equation (67) with the one-dimensional tight-binding Hamiltonian ˆ H := − X x ∈ Z (ˆ a † x ˆ a x +1 + ˆ a † x +1 ˆ a x ) (68) and the Lindblad op erator representing dephasing noise ˆ L x := p 2 γ dep ˆ a † x ˆ a x . (69) Under this setting, as shown in [42], the Liouvillian L can b e mapp ed to the Hubbard Hamiltonian with imaginary in teraction. T o sho w this, define sup erop erators, C x, ↓ [ ˆ ρ ] := ˆ a x ˆ ρ, (70) C x, ↑ [ ˆ ρ ] := h Y x ∈ Z e iπ ˆ a † x ˆ a x i ˆ ρ ˆ a † x , (71) and the Hilb ert-Schmidt inner pro duct for op erators ˆ A and ˆ B , ( ˆ A, ˆ B ) := T r[ ˆ A † ˆ B ] . (72) Then, the adjoint sup erop erators of C x, ↑ and C x, ↓ are given b y C † x, ↓ [ ˆ ρ ] = ˆ a † x ˆ ρ, (73) C † x, ↑ [ ˆ ρ ] = h Y x ∈ Z e iπ ˆ a † x ˆ a x i ˆ ρ ˆ a x . (74) These sup erop erators satisfy the canonical anti-comm utation relations, {C x,σ , C y ,σ ′ } = 0 , {C x,σ , C † y ,σ ′ } = δ x,y δ σ,σ ′ . (75) Therefore they can b e regarded as the creation and annihilation op erators of spin-1/2 fermions. The v acuum state annihilated b y all C x,σ corresp onds to the pro jector onto the original v acuum: ˆ Ω := | 0 ⟩⟨ 0 | . (76) 19 In terms of the ab o ve sup erop erators, the Liouvillian is equiv alent to the Hubbard Hamiltonian (8) with imaginary interaction: i U LU † = − X x ∈ Z ,σ ∈{↑ , ↓} ( C † x,σ C x +1 ,σ + C † x +1 ,σ C x,σ ) + 2 iγ dep X x ∈ Z (2 C † x , ↑ C x, ↑ C † x, ↓ C x, ↓ − C † x , ↑ C x, ↑ − C † x, ↓ C x, ↓ ) , (77) where we define the unitary sup erop erator as U := Y x ∈ Z e iπ C † 2 x − 1 , ↓ C 2 x − 1 , ↓ . (78) F rom this corresp ondence, the time evolution of a density matrix element is equiv alent to the propagator of the Hubbard mo del with imaginary interaction: ⟨ x 1 , · · · , x M | e L t [ | y M +1 , · · · , y N ⟩⟨ y 1 , · · · , y M | ] | x M +1 , · · · , x N ⟩ = h M Y j =1 ( − 1) x j − y j i e − 2 γ dep N ψ t ( x ; a 0 | y ; a 0 ) , (79) where a 0 := ( ↓ , · · · , ↓ | {z } M , ↑ , · · · , ↑ | {z } N − M ) . (80) Th us once the exact expression of the propagator is obtained, one can exactly calculate the time dep endence of the density matrix for finite-particles in a tigh t-binding c hain with dephasing noise. F urthermore, as shown in [42], correlation functions satisfy the same equation of motion as the densit y matrix elemen ts. Therefore, physical quantities for infinitely many particles can b e calculated from the propagator of finite particles. W e remark that a non-unitary circuit, which pro vides an in tegrable trotterization of this mo del, has b een prop osed in [53]. A.2 Hubbard mo del under tw o-b o dy loss The Hubbard mo del sub ject to t wo-bo dy loss is describ ed b y the GKSL equation (67) with the one-dimensional F ermi-Hubbard Hamiltonian (1) ˆ H = − X x ∈ Z X σ = ↓ , ↑  ˆ c † x,σ ˆ c x +1 ,σ + ˆ c † x +1 ,σ ˆ c x,σ  + 4 u X x ∈ Z ˆ n x, ↓ ˆ n x, ↑ (81) and the Lindblad op erators representing the tw o-b o dy loss ˆ L x := p 4 γ loss ˆ c x, ↑ ˆ c x, ↓ . (82) As sho wn in [43], the eigenv alues and eigenstates of this mo del are determined b y the Hubbard mo del with complex interaction, ˆ H eff := ˆ H − i X x ∈ Z ˆ L † x ˆ L x (83) = − X x ∈ Z X σ = ↓ , ↑  ˆ c † x,σ ˆ c x +1 ,σ + ˆ c † x +1 ,σ ˆ c x,σ  + 4( u − iγ loss ) X x ∈ Z ˆ n x, ↓ ˆ n x, ↑ . (84) 20 F urthermore, the probabilit y that N particles initially lo cated at sites y = ( y 1 , · · · , y N ) with spin configurations b = ( b 1 , · · · , b N ) are found at sites x = ( x 1 , · · · , x N ) with spin configurations a = ( a 1 , · · · , a N ) at time t is given by the propagator of the Hubbard mo del with the complex interaction (84), Pr[ x ; a | y ; b ] = | ψ t ( x ; a | y ; b ) | 2 . (85) This fact can b e derived from the decomp osition of the Liouvillian in to t wo parts, L = K + A [44, 45]: K [ ˆ ρ ] := − i ˆ H eff ˆ ρ + i ˆ ρ ˆ H † eff , A [ ˆ ρ ] := 8 γ loss X x ∈ Z ˆ c x, ↑ ˆ c x, ↓ ˆ ρ ˆ c † x, ↓ ˆ c † x, ↑ . (86) Note that K conserves the total num b er of particles, whereas A decreases it. B Inhomogeneous XXX structure of the scattering amplitude In this app endix, we describ e the relation b etw een the scattering amplitude (14) and the Bethe ansatz wa ve function of the inhomogeneous XXX spin chain, and deriv e (19) based on this relation. As explained in section 3, the scattering amplitude for charges given in (14) is the Bethe ansatz wa ve function of the inhomogeneous XXX spin chain [1, 46]. T o sho w this, w e first introduce the L -matrix at site n ( n = 1 , . . . , N ) as L n ( λ − µ ) := 1 λ − µ + iu  λ − µ + iuσ z n 2 iuσ − n 2 iuσ + n λ − µ − iuσ z n  , (87) where σ + n , σ − n , and σ z n denote the spin-1/2 raising, low ering, and z -comp onen t P auli matrices at site n , resp ectively . The mono dromy matrix of the inhomogeneous XXX mo del is then defined by T ( λ | s P ) := L N ( λ − s P ( N ) ) · · · L 1 ( λ − s P (1) ) . (88) Note that s j is giv en by (15), and it represents the inhomogeneity in the ab ov e equation. T o gain an intuitiv e understanding of the L -matrix and the mono dromy matrix, it is helpful to introduce their graphical representations. The L -matrix is depicted as the crossed lines shown in figure 5 (a). In this figure, the v ertical line with the rapidit y µ corresp onds to the quantum space, describing the spin-1/2 degree of freedom at site n whereas the horizontal line with the rapidit y λ represents the auxiliary space. The mono drom y matrix is represented by a single horizontal line in tersecting N v ertical lines, as illustrated in figure 5 (b), since the monodromy matrix is defined as an ordered pro duct of L -matrices. The B -op erator is given by the (1,2)-entry of the mono dromy matrix: B ( λ | s P ) := T ( λ | s P ) 1 , 2 . (89) Acting on the all-up state with these B -op erators yields the Bethe ansatz state of the inhomogeneous XXX mo del, B ( λ 1 | s P ) · · · B ( λ M | s P ) | ↑ , · · · , ↑⟩ . (90) 21 Figure 5: Graphical represen tations of the L -matrix and monodromy matrix. (a) L -matrix L n ( λ − µ ), (b) Mono dromy matrix T ( λ | z P ) for the case of N = 4. Then one sees that this state is equiv alent to (14): | s P ; λ ⟩ = Y 1 ≤ m