A Riemann-Hilbert formulation for the finite temperature Hubbard model

Inspired by recent results in the context of AdS/CFT integrability, we reconsider the Thermodynamic Bethe Ansatz equations describing the 1D fermionic Hubbard model at finite temperature. We prove that the infinite set of TBA equations are equivalent to a simple nonlinear Riemann-Hilbert problem for a finite number of unknown functions. The latter can be transformed into a set of three coupled nonlinear integral equations defined over a finite support, which can be easily solved numerically. We discuss the emergence of an exact Bethe Ansatz and the link between the TBA approach and the results by J\"uttner, Kl\"umper and Suzuki based on the Quantum Transfer Matrix method. We also comment on the analytic continuation mechanism leading to excited states and on the mirror equations describing the finite-size Hubbard model with twisted boundary conditions.


Introduction
The Hubbard model [1] arises as an approximate description of correlated electrons in narrow-band materials. Despite its simplicity, it is believed to capture important nonperturbative features of real many-body fermionic systems. In one dimension, the model is exactly solvable and its relevance for the study of strongly correlated electrons is perhaps comparable to that of the Ising model for magnetism. Quite unexpectedly, the Hubbard model has also attracted the attention of high energy physicists due to its multiple connections [2,3,4,5,6,7,8] with the integrable spin chains emerging in the context of N =4 Super Yang-Mills theory [9,10,11]. In particular, the exact Bethe Ansatz equations of Lieb and Wu [12], and the corresponding string hypothesis [13], have many features in common with the asymptotic Bethe Ansatz of Beisert and Staudacher for the anomalous dimensions of single-trace operators with large quantum numbers [11].
We start from the 1D Hubbard Hamiltonian written in the form where L is the length of the chain, c † i,σ and c i,σ are fermionic creation-annihilation operators satisfying with periodic boundary conditions c 1,σ = c L+1,σ , c † 1,σ = c † L+1,σ , and u ≥ 0 is a dimensionless coupling constant proportional to the electric charge. For L even, the Hamiltonian (1.1) is SO(4)-symmetric [14]. We shall consider a two-parameter deformation of (1.1) where this symmetry is explicitly broken by coupling the electrons to a chemical potential µ and to a magnetic field B: with the conserved electron numberN and spin S z operators defined aŝ Lieb and Wu showed that this system is integrable and that it can be solved by means of the Bethe Ansatz (BA) method [12]. The spectrum of the Hamiltonian (1. and the corresponding energy is given by (1.7) The system at finite temperature can be studied using the Thermodynamic Bethe Ansatz (TBA) method [15,13]. The TBA equations for the Hubbard model, an infinite set of coupled nonlinear integral equations, were first derived by Takahashi starting from equations (1.5), (1.6) and the classification of their solutions according to the so-called string hypothesis [13]. An alternative approach was introduced more recently by Jüttner, Klümper and Suzuki in [16], and is based on an associated integrable lattice model introduced by Shastry [17,18] and the path integral approach to thermodynamics, also called Quantum Transfer Matrix (QTM) method [19,20,21,22]. One of the main results of [16] is a much simpler set of only three coupled nonlinear integral equations (NLIEs) of Klümper-Batchelor-Pearce-Destri-DeVega type [23,24,25,26]. Although the two approaches lead to equivalent expressions for the Gibbs free energy, a direct link between the two sets of nonlinear integral equations was until now not found for the Hubbard model. For the discussion of the equivalence between QTM and TBA in other models, see [27].
In this paper we shall reconsider the Hubbard model TBA inspired by some notable recent results obtained studying the TBA equations for N =4 Super Yang-Mills (SYM). The TBA was introduced in the AdS/CFT context [28,29,30] to overcome the so-called wrapping problem [31] affecting the Beisert-Staudacher equations. Recently, this very complicated set of TBA equations was recast into the greatly simplified form of a nonlinear matrix Riemann-Hilbert problem: the Quantum Spectral Curve or Pµ-system [32,33]. This new formulation, contrary to the TBA, treats the full spectrum on an equal footing. It has already led to an impressive number of perturbative and nonperturbative results [34,35,36]. A similar reduction was also recently obtained [37] in the case of the N =6 Chern-Simons theory, and made possible the determination of the so-called slope and interpolating functions, both nontrivial nonperturbative quantities [38]. The role of the Quantum Spectral Curve approach in the general integrable model framework is not fully understood and one of the purposes of the present work is to investigate whether a similar structure arises also in the context of the Hubbard model.
(1. 11) In (1.8) and (1.9), P α i denotes the second-sheet evaluation of P α i . The Gibbs free energy f is contained in the large-z asymptotics, for instance ln T H 1,1 (z) ∼ − 1 T (f + µ + u). (1.12) We also found an alternative formulation, based on the fact that the P functions satisfy the functional relation where F ab (z) (a, b = +, −) are entire functions analytic on the whole complex plane. The set of equations (1.13) and (1.8), (1.9) is equivalent to (1.8)-(1.10), and may be seen as the analogue of the AdS/CFT Pµ-system. Furthermore, these relations imply that the zeros of the P and F functions are constrained by a set of exact Bethe Ansatz equations. More precisely, setting Q + (z) = e iBz/u P H + (z) P H + (z) and Q +− (z) = e i 2 (B−μ)z/u F +− (z), we have Notice that the solutions of these equations have an infinite number of Bethe roots. In fact, (1.14), (1.15) coincide with the infinite Trotter number limit of the exact Bethe Ansatz diagonalising the Quantum Transfer Matrix. Since the latter equations are the starting point for the derivation of the NLIEs of [16], our analysis provides the missing link between the two different approaches to the Hubbard model thermodynamics. As mentioned above, this was a longstanding problem (see, for example, the discussion at the end of Chapter 13 of [39]).
In this paper we also present a preliminary study of the numerical solution method for the ground state. At least for sufficiently high values of the temperature, we expect that the ground state is singled out by the requirement that all the roots s i are located on the second sheet. In the region of validity of this assumption, our functional relations give rise to a simple set of nonlinear integral equations, easy to solve numerically. We have explicitly checked the agreement of our results against the TBA predictions for many points in the region |B| < 1, |µ| < 1, 0 < u < 2 and T ≥ 1. However, the study of [16] indicates that there is a flow of roots from the second to the first sheet as T is decreased. It should be possible to generalise our numerical method to this parameter region, as well as to excited states, but we leave this problem for future studies.
As a last remark, we point out that this infinite Bethe Ansatz is reminiscent of the one discovered in the AdS/CFT context in [40,41] (however, in the latter case a relation with a lattice construction is not known).
The rest of the paper is organised as follows. Section 2 contains the TBA equations written in a slightly modified form that highlights some of the symmetries of the model and can be implemented numerically without range restrictions on the parameters µ and B. Section 3 contains the Y-system [42] and the basic discontinuity relations [43] that allow its extension to arbitrary branches of the Riemann surface on which the Y functions, solution of the TBA equations, are defined. In Section 4, taking hints from the recent results in AdS/CFT and motivated by the symmetry of the TBA equations, the Y's are parametrised in terms of T functions in two alternative gauges (T V and T H ). The simple analytic properties in these two gauges motivate a further reparametrisation of the T's as 2×2 determinants of more elementary objects: the P functions. Moreover, a resolvent-type representation for the P's allows to express all the relevant quantities in terms of a pair of densities ρ V and ρ H fulfilling a new set of NLIEs defined on a finite support (described in Section 7). Equation (1.13) and the exact Bethe Ansatz equations are derived in Section 5, where the connection with the results of [16] is also discussed. The free-fermion limit is briefly discussed in Section 6. Section 7 contains the new NLIEs together with preliminary numerics and a brief description of the numerical technique adopted. Section 8 describes the formal transformation relating the thermodynamic equations to the finite-size ones describing the Hubbard chain with twisted boundary conditions. The more technical parts of the analysis are confined to four Appendices. A proof of the special analytic properties of the Y functions on the Riemann section with only short cuts (the "magic" sheet of [44]) is reported in Appendix A. The study of the monodromy properties of the P functions and the proof that they live on a two-sheet Riemann surface is given in Appendix B. The derivation of the factors connecting the T V to the T H gauge, the proof of equation (1.10) and the derivation of the novel set of NLIEs are the main results of Appendix C. Finally, Appendix D contains a dictionary linking this work to the paper [16].

Thermodynamic Bethe Ansatz equations
In the notation of Takahashi, the solutions of the TBA equations are η n (z), η n (z), ζ(k) , n ∈ N + , (2.1) where the physical range for the arguments is z ∈ R and k ∈ [−π, π]. It is very convenient to reparametrise the variable k as z = sin(k), and introduce the two-indexed functions Y m,n as follows: The structure of the TBA equations is thus codified on the L-shaped diagram represented in Figure 1, with every node of the diagram associated to one of the unknown Y functions of the TBA. As a consequence of the change of variable z = sin(k), the functions Y 1,1 (z) and 1/Y 2,2 (z) have a branch cut of square root type on the real axis for z ∈ (−1, 1), and are in fact two branches of the same function ζ(z): denoting with a tilde the analytic continuation around one of the branch points z = ±1, we have Y 1,1 (z) = 1/Y 2,2 (z).  The TBA equations consistent with the Hamiltonian defined in (1.3) can be written as where we have defined In (2.5)-(2.7), the convention for the convolutions is a * b(z) = R dv a(v)b(v − z). Furthermore, here and in the rest of the paper we are implicitly assuming that the integrals over the interval z ∈ (−1, 1) run slightly above the branch cut. Notice that the TBA equations (2.5)-(2.7) are written in a slightly different form as compared to [13]. For the equations of [13], the standard method of iterative solution [45] requires sign restrictions on B and µ, while here this can be avoided since the symmetries B ↔ −B and µ ↔ −µ of the ground state solution are explicit. A further symmetry corresponds to the exchange of the two wings up to a change of sign for all the driving terms (T ↔ −T, φ ↔ −φ).
There are many equivalent expressions for the Gibbs free energy, for example: The free energy can equivalently be rewritten using only the Y functions of the horizontal part of the diagram of Figure 1, as (2.14) The equivalence between (2.12) and (2.14) can be checked by using the TBA equation (2.7), and reflects the symmetry (2.11) and the properties [39]: Equations describing excited branches of the free energy -which control the correlation lengths among local operators at finite temperature -can in principle be obtained by analytic continuation [46], see also [47,48,49,50]. For the Hubbard model, this has been accomplished only for a few states [51,52,53]. The TBA equations appear not to be the optimal tool for this analysis, and the results of [51,52,53] were obtained adopting the QTM method. It would be important to have a more complete understanding of the free energy spectrum, and one of the aims of the present work is to propose an alternative formulation which seems to have some advantages over the existing approaches.

The Y-system and the discontinuity relations
In this Section we shall describe the essential analytic properties of the system, necessary for the simplifications to be . An important complication comes from the fact that the Y functions have square root branch points at certain positions in the complex rapidity plane.
In particular, it can be seen analysing the TBA equations that Y 1,n and Y n,1 are analytic in a strip |Im(z)| < (n − 1)u, but have square root branch points at z = ±1 + i(n − 1)u, z = ±1 − i(n − 1)u, and possibly at other positions further from the real axis.
To describe this structure, let us setup some useful notation. We adopt the convention that, for complex z, the Y's are evaluated on a Riemann sheet with short cuts, the "magic" sheet (see Figure 3), and use the notation g(z) for the analytic continuation of a function g(z) around either 1 of the branch points z = ±1. Since we reserve this notation to branch points of square root type, we always have the property g = g. Another notation that we shall use often is g [+n] (z) = g(z + inu), where the shifts are evaluated on the magic sheet.
From (2.5)-(2.7) with n = 1, one can derive the following discontinuity relations [43,54] for the monodromies around some of the branch points closest to the real axis Furthermore, a crucial property implied by the TBA equations is that the Y's are solutions 1 It can be checked that the result of the analytic continuation around two branch points symmetric with respect to the imaginary axis is always the same.
of the Y-system [42] Y m,n (z + iu)Y m,n (z − iu) = (1 +Y m,n+1 (z))(1 +Y m,n−1 (z)) (1 + 1/Y m+1,n (z))(1 + 1/Y m−1,n (z)) , (3.4) where (m, n) ∈ {(1, k), k ∈ N + }∪{(k, 1), k ∈ N + } with boundary conditions Y k,0 = 1/Y 0,k = Y k+2,2 = 1/Y 2,k+2 = 0 for k ∈ N + (see Figure 1). Notice that there is no independent equation centered at the node (2, 2). Indeed, the function Y 2,2 is simply related to Y 1,1 through equations (3.1). In (3.4), we have used the notationY to emphasise that the functional relations (3.4) are valid on the specific Riemann section shown in Figure 4, the "mirror" sheet, where all the branch cuts are traced as semi-infinite lines parallel to the real axis, and do not cross the strip |Re(z)| < 1. This gives a precise prescription on how to evaluate the shifted values appearing on the lhs of (3.4). Adapting the arguments of [43,54,55], it can be proved that the functional relations (3.1)-(3.4), together with the asymptotics 6) and the assumption that the Y's have no zeros or poles in the strip |Im(z)| ≤ u, are fully equivalent to the TBA. One might expect that each Y function should display, on a generic Riemann section, an infinite ladder of further square root branch points at steps of 2iu, replicated from the ones closer to the real axis by the Y-system (3.4) [43]. However, this is not the case. We shall indeed prove that each of the Y functions has only a finite number of branch points on any sheet and that the number of Riemann sheets is actually finite. The fact that all functions appearing in this problem are defined on a finite genus Riemann surface is perhaps obvious from the perspective of the exact Bethe Ansatz and the Quantum Transfer Matrix construction of [16]. However, starting from the TBA equations this property is much harder to prove. To establish this result, we adopt the following strategy. First, we show that the Y functions have at most four branch cuts on the magic sheet: • Y 1,n (z) and Y n,1 (z) for n ≥ 2 have only four branch cuts at z ∈ (−1, 1) ± iu(n − 1), z ∈ (−1, 1) ± iu(n + 1), • Y 1,1 (z) and Y 2,2 (z) have only three branch cuts at This result is established in Appendix A 2 . Secondly, it is shown in Appendix B that no further branch cuts can appear on the other sheets. A direct numerical solution of the TBA 2 It is important to notice that the proposed proof can be straightforwardly adapted to the known AdS/CFT cases, this gives an alternative and perhaps more transparent way to understand the nice properties discovered in [44] for the AdS5/CFT4 Y functions on the magic sheet.
are displayed in Figures 5 and 6.

The T-system
Following the strategy which allowed the simplification of the AdS/CFT TBA, we shall now make a chain of simplifications and reduce relations (3.1)-(3.4) to the finite set of constraints (1.8)-(1.10). The first step is to parametrise the Y functions as where the T's (T n,s with n, s ∈ N) are defined on the extended lattice displayed in Figure  2 and satisfy, on the mirror section, a discrete Hirota equation (the T-system) is automatically fulfilled as a consequence of (4.1) and (4.2). However, the parametrisation (4.1) is not one-to-one. Different solutions of the T-system, corresponding to the same set of Y's, are linked by a gauge transformation [56]. A clever gauge choice can greatly simplify the problem.

Vertical and horizontal gauges
We shall now introduce the two gauges T H and T V , each defined by a set of simple conditions on one of the wings, horizontal (H) or vertical (V), of the diagram in Figure 2. We require that they fulfill the following properties, compatible with the cut structure of the Y functions, summarised in Section 3: T H 1,n (z), (n ≥ 1) has only two branch cuts at z ∈ (−1, 1) ± inu (4.3) and T V n,1 (z), (n ≥ 1) has only two branch cuts at z ∈ (−1, 1) ± inu (4.5) The conditions listed above still leave a residual gauge invariance. These gauges are fixed uniquely by imposing the following extra conditions: It is shown in Appendix B that this choice is indeed possible. Given these definitions, it is far from obvious that also the functions T H n,m with n > m and T V n,m with m > n have such a simple cut structure. However, one of our main results, is that these two gauges are in fact connected by an elementary transformation: To prove (4.8)-(4.10), it is necessary to parametrise (4.3)-(4.6) with more elementary building blocks, the P functions introduced in Section 4.2 and the resolvent parametrisation of Section 4.3. The technical details of the derivation are described in Appendix C.
In the matching condition (4.8)-(4.10), together with the requirement that the T's have no poles, is hidden the full content of the original TBA equations, and its generalisation to excited states. Furthermore, we will see in Section 5 that the T functions have a clear interpretation: they are directly related to the eigenvalues of the (fused) Quantum Transfer Matrix in the infinite Trotter number limit [16].

The P functions
The next natural step is to parametrise the T functions defined in (4.3)-(4.7), and possessing only two cuts, in terms of objects having only one branch cut running along (−1, 1). We introduce a pair of P a (a = +, −) functions for each of the wings, and write with s ∈ N + . Thanks to relations of Plücker's type among determinants, this parametrisation is automatically consistent with the T-system in the corresponding H or V wing of the diagram. This "quantum Wronskian" construction [57,56] suggests that the P's are related to the Q functions describing the spectrum of the thermodynamic problem [16]; this relation is made precise in Section 5.
To fix completely these functions, we have to specify their asymptotics and to constrain further their analytic properties. As suggested by the asymptotics of the Y functions (3.5), we demand that 3 , for large z, The particular normalization of the prefactors C H and C V in (4.13) is chosen for future convenience. As shown in Appendix B, the four P functions live on a two-sheet Riemann surface. This is an important simplification as compared to the AdS/CFT case studied in [32], where the analogous P functions have additional infinitely many branch cuts, starting from their second sheet 4 . As we discuss in Appendix C, equation (4.7) gives the constraint Equations (4.14)-(4.15), together with the relation giving the matching of the two wings: the asymptotics (4.13) and information on the number of zeros on the first Riemann sheet, form a closed set of conditions for the P functions. In Section 7, we will show how to solve these equations numerically. Notice that the constants C H and C V are fixed by the solution to this system, and contain the free energy f through the relation: We will prove (4.17) in Section 4.4, with the aid of a very convenient parametrisation in terms of resolvents, inspired by [44,41].

Resolvent parametrisation
In the following, we shall assume that the P functions have no zeros on the first sheet. We expect that this singles out the ground state solution in a large parameter region for sufficiently high temperatures. The numerical solution of the functional relations presented in Section 7 reveals the existence of an infinite number of zeros on the second sheet.
Starting from the ground state, we expect that excited states can be obtained by analytic continuation in B and µ, leading to the migration of a finite number of zeros to the first sheet and to simple modifications in the equations presented in this Section. We point out that, in a neighbourhood of T = 0, even the ground state solution should present a finite number of zeros on the first sheet [16], and therefore the following analysis needs to be modified.
Let us introduce two resolvent functions where the densities ρ H and ρ V are for the moment undetermined. This is the most generic parametrisation of a function with a single cut, no zeros and no poles on the first sheet, and large-z asymptotics G α ∼ 1/z on the first sheet. Under analytic continuation through the cut, the resolvents transform as Thus, under the assumptions discussed above, we can parametrise the P's as where the functions h α (α = H, V) are required to have a single cut, no zeros or poles and a constant leading asymptotics on the first sheet. Taking ratios of these functions we find As a consequence of the simple monodromy properties of the P's, also ρ H and ρ V live on a two-sheet Riemann surface, with More precisely, we require that the behaviour at the branch points is such that the functions ρ α (z)/(z 1 − 1/z 2 ) are analytic in a neighbourhood of the cut z ∈ (−1, 1). Due to (4.21), however, the densities can have logarithmic singularities elsewhere in the complex plane, in correspondence to zeros of the P's. Substituting (4.20) into (4.14) and (4.15), we find  The solution of the latter equations, compatible with the asymptotics (4.13) and the absence of zeros on the first sheet, is unique and given by 5 : Thus, with the definitions (4.25) and (4.26), the parametrisation (4.20) automatically fulfills the constraints (4.14),(4.15).
Considering the large-z asymptotics of (4.25),(4.26) we find that, on the first sheet, In Section 4.4, we will connect the integrals appearing in (4.27) to the free energy and establish relation (4.17). From (4.11)-(4.12) and (4.20), one can rewrite the T functions in terms of the densities as where we have borrowed the notation of [44,41] and denoted (4.30) 5 Notice that in (4.25), φ(v) needs to be evaluated just above the branch cut, so it agrees with 1 Equation (4.29) defines a gauge transformation between the T α and the T α functions, with In particular, the h factors in (4.29) cancel in the gauge invariant combinations (4.1). The gauges T α are characterised by the simple large z asymptotics which fixes them uniquely. Finally, in these gauges relations (4.7) become As discussed in Section 7, the densities can be computed by solving numerically a set of NLIEs, which can be derived from the constraint (4. 16). Some examples of solutions are displayed in Figures 7 and 8. Finally, let us remark again that the information contained in the pair {ρ H , ρ V } is fully equivalent to the knowledge of the solution of the TBA equations. In Figure 9, the TBA solution is compared with the same quantity reconstructed from the parametrisation (4.1),(4.30).

The free energy
The purpose of this Section is to derive a simple formula for the free energy in terms of the densities. The starting point is equation (2.12), written in the equivalent form: where we have setT H 1,s = T H 1,s sinh(B)/(sinh(sB)), so that lnT H 1,s (z) ∼ 0 at large z. Adopting the "telescoping" technique used in [55,44,40], it is now possible to remove the infinite sum appearing on the rhs of (4.35). SinceT H 1,n is analytic for |Im(z)| < (n + 1)u, we can shift the integration contours and prove that, for n ∈ N + , This shows that many cancellations take place in (4.35), leading to Using (4.1) and the fact that Y 2,2 = 1/ Y 1,1 , (4.37) becomes The convolution appearing in the first line of (4.38) can be viewed as a contour integral around the cut of k (z), and, deforming this contour to a pair of infinite lines, one can write Shifting the line contours appearing in the latter expression (and pushing some of them to infinity), one can show that the integrals involving T H 1,1 in equation (4.38) actually cancel each other out. Substituting (4.34) for T H 1,0 , we finally arrive at the desired result: Starting from the alternative expression (2.14) for the free energy, working in the gauge T V and following the same steps, we also obtained The comparison between (4.40), (4.41) and (4.27) finally gives the formula (4.17), showing that the free energy f appears in the asymptotics of the P's.

Energy-carrying Bethe roots
Let us look more closely at the kernels appearing in (4.40) and (4.41). From (4.21), we have where we have used (4.14). The same relation withB →μ, H → V is valid for the density ρ V . Thus, the kernels entering the energy formula have logarithmic singularities in correspondence to the zeros of the P functions. Should any of these zeros cross the integration contours in (4.40),(4.41), the energy would pick extra residue terms of the form iT k(z i ), as is easy to verify by first performing an integration by parts 6 . This suggests that zeros of the P functions can be regarded as energy-carrying Bethe roots for the finite temperature spectrum. According to our assumptions, for the ground state with real B and µ all the zeros lie on the second sheet, but they may cross the interval and move to the first sheet under analytic continuation in B and µ. We will return on this point in Section 6 below.

Exact Bethe Ansatz and relation with the Quantum Transfer Matrix
In this Section we will recast the closed set of conditions (4.14)-(4.16) into another, interesting form, which reveals the presence of an exact Bethe Ansatz. This second formulation is perhaps closer in spirit to the Pµ-system of AdS/CFT [32], which was one of the main sources of inspiration for the present work.

Formulation as a coupled Riemann-Hilbert problem
We start by considering the analytic continuation of (4.16) through the cuts at z ∈ (−1, 1)± iu, which yields Two of the P functions in (5.1) can be eliminated using the constraints (4.14), (4.15). Solving (5.1) for the remaining two P's, we find where F ab has a simple explicit expression: Relation (5.3) shows that F ab (z) could -in principle -have a pair of branch cuts on the lines Im(z) = ±iu. However, by shifting equation (5.2) of ±iu, and evaluating the discontinuity across the branch cut, we find which implies that F ab is analytic in the whole complex plane. The absence of branch cuts can also be directly verified by using the definition (5.3): the quantities F ab vanish since they are proportional to the difference of the lhs and rhs of (4. 16). Notice that the lhs of (5.2) has no poles. As a consequence, it is possible to prove that the F functions, appearing on the rhs, are entire on the whole complex plane. This can be also deduced from equation (5.9) below.
These newly introduced quantities fulfill, together with the P's, a closed set of functional relations. The fundamental set of equations is with the requirement that the F functions have no cuts and no poles, and the P's have only one cut on each of their two sheets, and no poles. These equations make explicit the nonlinear Riemann-Hilbert nature of the problem, since they show how the two branches of the P's are connected through the entire 2 × 2 matrix F ab . In this respect, they are reminiscent of the Pµ-system. It is important to underline that (5.5)-(5.7) are completely equivalent to (4.14)-(4.16). In particular, it is possible to show that the definition (5.3) and the relation (4.16) are both hidden in these equations. Equations (5.5)-(5.7) imply the existence of many other functional relations among the F's and the P's. For instance: Finally, as already remarked in Section 1, the Hubbard Hamiltonian (1.1) has, for B = µ = 0, a hidden SO(4) ∼ = SU (2) × SU (2)/Z 2 symmetry [14]. While this symmetry was not very evident in the original Y-system and discontinuity relations, it appears to be nicely encoded in the structure of equations (5.5)-(5.7). At generic values of the magnetic field or chemical potential, the symmetry is broken by the boundary conditions (4.13).

The exact Bethe Ansatz
The relations obtained in the previous Section contain an exact Bethe Ansatz constraining the position of the zeros of the P and F functions. Starting from (5.2), we immediately get , a, b = +, −. (5.9) Combining these expressions to form bilinear combinations of the F ab 's reveals a further set of interesting relations of quantum Wronskian type. For instance, using (5.9) and (4.14),(4.15), we obtain The rhs of (5.10) defines an entire function on the whole complex plane. It is convenient to isolate the exponential prefactors and set The a = +, b = − case of equation (5.2) can then be rewritten as We assume that for the ground state all zeros of the P functions are on the second sheet; therefore the zeros of the rhs of (5.12) coincide with the zeros of Q H + , giving the first of the Bethe Ansatz equations quoted in Section 1. The second BA equation can be derived by taking ratios of the expressions obtained by evaluating (5.12) at zeros of Q , at Q +− (w α ) = 0. (5.14) Very interestingly, these equations appear to be the infinite Trotter number limit of the BA diagonalising the Quantum Transfer Matrix of Jüttner, Klümper and Suzuki [16]. We remind the reader that the QTM is defined as a discrete object acting on aN -site quantum space (N is known as the Trotter number), and that the thermodynamics of the Hubbard model is described by the largest eigenvalue of the QTM in the limitN → ∞. For the ground state at finite even values ofN , the BA of [16] is where 7 limN →∞ bN (z) = e −2φ(z) and The functions appearing in the BA equations (5.13),(5.14) are the continuum version of these polynomials. In fact, we believe that they can be factorised over their zeros as infinite products of the form 8 where Re(s −i ) = −Re(s i ), Re(w −α ) = −Re(w α ). In (5.19), the zeros have been paired up in order to make the product convergent. This rearrangement is necessary, since the zeros accumulate at infinity with an evenly spaced asymptotic distribution (see Sections 6 and 7). The regularisation (5.19), without extra Hadamard factors, is consistent with the infinite Trotter number limit of (5.15)-(5.17). In addition to (5.13), (5.14), there are naturally other equivalent sets of BA equations. For completeness, let us discuss the general structure. From (5.9) and (4.14), (4.15), one can obtain the following quantum Wronskian-type relations whereB ± = ±B,μ ± = ±μ. We expect that all the Q ab and Q a functions thus introduced admit a factorisation of the form (5.19), each with a different set of zeros. The systems of 7 See Appendix D for more details. 8 We expect that the P functions can be factorised in the following form (cf [41]), where z  Finally, let us point out that the definitions (5.24), (5.25) can be inverted as (a = +, −), leading to the following equivalent formulae for the free energy: (5.32)

Relation with the Quantum Transfer Matrix
By comparison with the results of [16] we can uncover the physical meaning of T H 1,1 , T V 1,1 and show that these functions are simply related to the eigenvalues of the Quantum Transfer Matrix. Starting from (5.2), we get , a, b = +, −.
The precise details of this identification are given in Appendix D. Notice that, as a consequence of the absence of poles for the P functions, this quantity does not have poles on any sheet, and that this pole-free condition gives precisely the Bethe Ansatz. Furthermore, setting P H a P H a =Q H a , we can write In the quantity in the square brackets we recognise the infinite Trotter number limit of the "auxiliary" transfer matrix eigenvalues Λ aux introduced in [16], equation (25). Setting a = +, the identification is (5.39)

The free fermion limit
The exact solution of the TBA equations at u = 0 was found already by Takahashi in [13]. It is interesting to recover this result starting from relations (4.14)-(4.16). In the limit u → 0 + , the shifts in the P functions shrink to zero and, for Re(z) ∈ (−1, 1), the P functions collapse to their values above/below the cut. Therefore, and equation (4.16) becomes A second constraint is simply obtained by continuing (6.3) to the second branch: Solving (6.3),(6.4) with the aid of (4.14),(4.15), we find We can now compute the densities using (4.21). For the horizontal wing the result is and, denoting the rhs of this equality as e 2Bρ H (z) = R(B,μ; z), the density characterising the vertical wing is given by e −2μρ V (z) = R(μ,B; z). The kernel appearing in the free energy formula (4.40) reduces to ln sinh 2 (2φ(z)) sinh 2 (B) with σ 1 , σ 2 = ±1, giving the well-known result for the Gibbs free energy at u = 0. The pattern of Bethe roots displayed by the free fermion solution is interesting, as the numerical solution for u > 0 (see Section 7) suggests that the zeros are smoothly deformed away from their positions at u = 0. Each of the P functions has two infinite strings of zeros on the second sheet, corresponding to the two factors on the rhs of (6.5). Denoting as z n (B,μ) the solution of the equation 2 φ(z n (B,μ)) +B +μ = i(2n + 1)π, (6.10) and defining A(B,μ) = z n (B,μ), n ∈ Z , (6.11) the distribution of zeros is summarised in Table 1, and illustrated in Figures 10 and 11 for two of the P's at B = 1 and µ = 1/2. Zeros belonging to A(B,μ) accumulate at infinity along the line Im(z) = −B/2 − µ/2, and their asymptotic spacing is πT . As can be obtained from equations (5.20)-(5.23) in the u → 0 limit, in the free fermion case the zeros of the F functions are a subset of the zeros of the P's (see Table 1). Let us make a short comment on the analytic continuation mechanism governing the transition to excited states. As already anticipated in Section 4.5, the free energy acquires an extra residue whenever one of the energy-carrying Bethe roots crosses the integration contour. The movement of these zeros can be driven by analytic continuation in B and µ. In general, such a crossing does not correspond to a branch point in the domain of the parametersB orμ, since the contour can be deformed to avoid the contact with the Table 1: Distribution of the Bethe roots among different P and F functions, with a, b = +, −, andB ± = ±B,μ ± = ±µ. Zeros of P's live on the second sheet.
wandering zero. Genuine branch points correspond to the so-called pinching phenomenon, when a pair of zeros collide on the contour from opposite sides [46]. In the free fermion case, this happens at values of B and µ given by: where λ B = e B/T , λ µ = e µ/T are the fugacities, corresponding to a pair of zeros of sinh 2 (B)/ sinh 2 (Bρ H ) = sinh 2 (μ)/ sinh 2 (μρ V ) pinching the contour of integration at the origin in the z-plane. Analytic continuation around one of these points causes the transition to an excited branch of the free energy, and of the P functions. Finally, it is interesting to notice that the branch points (6.12) mark the boundaries of the four phases of the system at T = 0. In the interacting regime u > 0 at T = 0, the phase diagram includes a fifth phase describing the Mott insulator behaviour [12] (see also Chapter 6 of [39]). It would be very interesting to investigate the branching structure of the free energy at finite temperatures and coupling, and link it to the phase diagram of the Hubbard model. We expect this to be possible with the numerical method described in Section 7 and plan to come back to these questions in the future.

Numerical solution
In Section 4.3, we have parametrised the P functions appearing in the problem in terms of the densities ρ H and ρ V . They can be computed by solving a system of coupled nonlinear integral equations, which determine simultaneously the densities and the function Y 1,1 entering the TBA equations. This formulation is very similar to the set of NLIEs proposed in [44] for the AdS 5 /CFT 4 spectral problem. We discuss here only the ground state equations.

Nonlinear integral equations
First, by expressing Y 1,1 in the T H and in the T V gauge, we find where 2) and the functions T i 1,1 , i = H, V depend on the densities through the parametrisation (4.30). In an iterative scheme, equations (7.1) can be used to update the values of the two densities starting from the knowledge of r(z). The numerical method is described in [44] and is reviewed below in Section 7.2. To close the system, there is a further equation determining Y 1,1 as a function of ρ H and ρ V : This relation is derived in Appendix C, both from the TBA equations and also purely from the Riemann-Hilbert formulation (4.14)-(4. 16). An equivalent form of (7.3), which is convenient for the numerical implementation, is where − denotes Cauchy's principal value integral.

The numerical method
The system (7.1),(7.4) can be solved iteratively for the values of ρ H (z), ρ V (z) and Y 1,1 (z) on the interval z ∈ (−1, 1). One iteration step, updating the values of the densities ρ (7.5) After each step, we update the solution as ρ . The introduction of the weights θ, 1 − θ is a common recipe used to ensure convergence. In all cases we considered, the scheme was stable taking θ = 1 2 . In the above described procedure, the most difficult step is the solution of (7.1) for the densities. Let us review the method discussed in [44], concentrating on the equation for the horizontal wing. Using (4.30) and the basic property G H = G H + ρ H , the first equality in (7.1) can be written as where / G H denotes the Cauchy principal value integral By extracting 9 ρ H from the rhs of (7.6), we obtain the density in terms of G H , / G H and r. This equality is used to update the value of ρ H in the last passage of (7.5).
The numerical evaluation of the singular integrals appearing in the NLIEs (7.1),(7.4) can be performed very efficiently using a Chebyshev expansion. To optimize the numerical method we found it convenient to discretise the densities by using a Chebyshev expansion of the second kind (where we have taken the correct parity into account) and evaluate principal value integrals using the properties where T n and U n denote the Chebyshev polynomials of the first and second kind, respectively. To produce the data presented in this paper, we took N trunc = 50. In the vast majority of the cases we considered, less than 30 iterations were sufficient to achieve convergence of the coefficients c (α) n entering (7.8) on the fifth digit. We observe that the error is approximately halved at every iteration.
Let us also make a comment on the region of convergence. We observe that, for fixed values ofB,μ and u, the iterative scheme is convergent for sufficiently high temperatures -in particular a preliminary study suggests that, for arbitrary B, µ ∈ R, the convergence region probably includes 0 < u < 2 and T ≥ 1 -but breaks down below a certain threshold temperature. Lowering u, the breakdown temperature decreases, and this hints that, for a given value of T , the method should be applicable without modifications in a nonvanishing neighbourhood of the free fermion point. We suspect that the breakdown of the method for low temperature or strong coupling is related to the appearance of zeros on the first sheet for the ground state solution [16]. We plan to come back to this issue in the near future.

Exploring the complex plane
The numerical method we have described computes ρ H (z), ρ V (z) and Y 1,1 (z) for −1 < z < 1. Once we have a solution on the interval, we can reconstruct the behaviour of these functions in the complex plane. We do this in two steps. First, we use equation (7.3) to compute Y 1,1 (z) for z on the whole first Riemann section -with cuts at z ∈ (−1, 1), z ∈ (−1, 1)±2iu -from the values of the densities on the interval. From the same information, we can also compute G H (z) and G V (z) for an arbitrary complex value of z. Then, one can obtain the values of ρ H (z), ρ V (z) for any complex values of z by inverting (7.1). The complex zeros of sinh(ρ H (z)) are visible in Figure 12 for the numerical solution corresponding to µ = 1/2, B = 1, u = 0.1, showing a clear qualitative similarity with the free fermion case.

Mirror equations
In this Section we shall describe the finite-size versions of the functional relations (5.5)-(5.7), which form a set of equations fully equivalent to the Lieb-Wu quantisation conditions (1.5),(1.6) with more general twisted boundary conditions. Comparing the thermal BA (1.14),(1.15) and the Lieb-Wu equations (1.5),(1.6), and considering the dispersion relations implied by (1.7), (4.40), one can see that the two systems are formally related by a simple map swapping L ↔ 1/T and single-particle energies with momenta 10 : Applying (8.1) to (5.5)-(5.7), we find a simple set of functional relations: where we require that the functions p α a (z) live on a two-sheeted Riemann surface, while the f ab 's are entire. Relations (8.2)-(8.4) can be seen as a set of Baxter-like equations for the Hubbard Hamiltonian. For L an even integer, they admit many solutions where the f 's are polynomials; consequently, the p's can be written as polynomial functions of x(z) and x(z), fixed in terms of their zeros on two sheets. For completeness, we can also consider exponential prefactors of the same type as the ones in (5.11), relabeling B → α, µ → β. We then find easily the quantisation conditions for the following Bethe parameters: the zeros s j of the p's, and λ l of the f 's. Setting k(s j ) = k j , and following the same route of Section 5, we find which is precisely the BA diagonalising the Hubbard Hamiltonian on the L-site chain with twisted boundary conditions [58,59]: The relation between the two Bethe Ansatz systems discussed above is directly connected to the path integral approach to the thermodynamics. Adopting the notation of [60], the latter is based on rewriting the partition function, of a spin chain of L sites with Hamiltonian H, as where Z 2D classical (N , L, u) is the partition function of an appropriately defined two-dimensional statistical model living on a L ×N lattice. Its column-to-column transfer matrix is the Quantum Transfer Matrix T QTM (u). By a π/2 rotation on the lattice we switch from (8.8) to the description The exact Bethe Ansatz (1.14),(1.15) describes the spectrum of T QTM in the Trotter limit N → ∞. It is an interesting open problem whether one can define the QTM directly in this limit, and give it a meaningful physical interpretation as a continuum model living on a space of size 1/T . A somehow similar problem has arisen in the context of AdS/CFT integrability, where Zamolodchikov's ideas on the TBA for Lorentz invariant scattering theories [61] were adapted to the study of the non-relativistic string sigma model dual to planar N =4 SYM by considering its doubly Wick rotated counterpart (the mirror model ) [31,62]. Thanks to the knowledge of the action for the AdS 5 × S 5 string sigma model, it has recently been possible to identify the corresponding mirror model as a string theory living on a mirror background [63].
Finally, we would like to mention that is also possible to introduce a mirror version of Takahashi's TBA, in such a way that it is equivalent to the finite-size BA (8.5), (8.6). Apart from minor subtleties 11 , the mirror TBA equations can be obtained from (2.5)-(2.7) and the energy formula (2.12) through the following formal map : (8.11) and by modifying the integration contours as (cf [62,28,29,30]) so that the roles of the mirror and magic sheets of Figures 3, 4 are interchanged 12 . From these equations, one could in principle repeat the reduction presented in this paper and obtain the system (8.2)-(8.4).

Conclusions
The one-dimensional fermionic Hubbard model is one of the most interesting systems of lowdimensional condensed matter physics. Since its appearance in 1963, it has been intensively studied by means of exact and perturbative methods, greatly advancing the understanding of the physics of electron transport in 1D solids. A partial grasp about the huge number of results on this model can be obtained by consulting [65], the book [39] and the collections of works in [66,67]. The purpose of this article is to add a little piece to the jigsaw, by recasting the Thermodynamic Bethe Ansatz equations of Takahashi as a nonlinear Riemann-Hilbert problem, reminiscent of the Quantum Spectral Curve formulation recently obtained for the study of anomalous dimensions in AdS/CFT [32]. One of the main results presented in this paper is a new set of nonlinear integral equations describing the thermodynamics of the system. In their region of validity (see discussion at the end of Section 7.2), these equations can be integrated numerically with very high precision and, even when implemented on Mathematica, the iterative algorithm converges in only a few seconds of CPU time. The complexity of this formulation is comparable to the system of nonlinear integral equations derived by Jüttner, Klümper and Suzuki in [16]. However, as a consequence of the fact that the equations proposed here are defined on a finite support, we think that they may prove more convenient for the study of finite temperature correlation lengths.

A. The magic sheet
In this Section we will show that the Y functions have the following cut structure on the magic sheet: • Y 1,n (z) and Y n,1 (z) for n ≥ 2 have only four branch cuts: • Y 1,1 (z) and Y 2,2 (z) have only three branch cuts: z ∈ (−1, 1), z ∈ (−1, 1) ± 2iu.
As we will argue in Appendix B, a stronger property is true, namely the Y functions have no branch points outside the positions specified above, on any Riemann sheet. However, here we will present only the proof on the magic sheet, which is easier as it relies only on the structure of the Y-system and discontinuity relations. Before we start, we need to rewrite the Y-system (3.4) -originally defined on the mirror section -on the magic sheet. Using (3.1)-(3.3), it is simple to prove that the magic-sheet version of the Y-system is We will now test the cut structure of the Y 1,n+1 functions for n ≥ 1. First, let us introduce Notice that the term in brackets in (A.5) has no cut because it falls into the analyticity strip. Therefore, the cut structure of the lhs of (A.5) depends on the X n factors. We shall now show that X n have no cut on the real axis, for all n ≥ 2. Starting from n = 2, we can compute .
This expression manifestly does not have branch points on the real axis, as all terms on the rhs fall into their respective analyticity strips. A very similar calculation shows that and, literally by repeating this calculation with shifted indices, one finds the general case: Taking into account the analyticity strips of the functions on the rhs, this equation shows by induction that all X n are free of cuts. Therefore, the lhs of (A.5) is analytic in a neighbourhood of the real axis, meaning that Y 1,n with n ≥ 2 has possibly two cuts in the upper half plane at Im(z) = u(n ± 1), but no cut at Im(z) = u(n + 3). Moreover, because of the Y-system, no further cuts are possible and the Y 1,n 's with n ≥ 2 have only two short branch cuts in the upper half plane. By symmetry, the argument can be repeated for the lower half plane and for the Y n,1 functions with n ≥ 2. where R manifestly has no cut on the real axis and is defined by 1,2 ) 1 + (Y [2] 3,1 ) −1 /(1 + 1/Y [1] 2,1 ) (A.11) Because Y 1,1 (z)Y 2,2 (z) = e −4 φ(z) has no cuts outside the real axis, the only discontinuity can come from the factor: Using the discontinuity relations (3.1)-(3.3), we find This shows that Y 1,1 (and therefore also Y 2,2 ) does not have a branch cut with Im(z) = 4u. The Y-system, together with the results already obtained for the other Y functions, imply that no further cuts are possible.

B. Monodromy properties of the P functions
The purpose of this Appendix is to provide a proof for some statements made in Section 4.
1) First, we prove that it is possible to choose the T gauges in such a way that equation In the proof we will use the resolvent parametrisation of Section 4.3, which, rigorously speaking, is valid only for the ground state. However, we expect that the result holds in general.
3) Finally, we discuss how to infer that the second-sheet evaluation of the P functions, P, do not have other branch cuts apart from z ∈ (−1, 1). This shows that the P's live on a Riemann surface with only two sheets.
To prove these properties, we will adapt many of the arguments of [44] to the present case.
Proof of equation ( .

(B.5)
The equation involving Y 1,2 is automatically satisfied by the parametrisation (B.2), while the equation for Y 2,1 will be used later. Now, we compare the T-system equations at the nodes (1, 1) and (2,2). On the magic sheet, these two equations read The set of equations (B.17) can be used to prove that the P's live on a two-sheeted Riemann surface. By definition, these functions have a single cut on their first sheet, and parametrise some of the T H 's as Obviously, by symmetry all these results are valid both in the upper and in the lower half complex plane, and show that P H ± have no cuts outside the real axis.
C. Fixing the gauge factor T H /T V The main purpose of this Appendix is to derive the form of the gauge transformation (4.8)-(4.10) relating the T H and T V functions. Along the way, we will also obtain the two NLIEs (C.4) and (C.8), adopted for the numerical solution method described in Section 7. Let us briefly discuss the logic of this Section. In Appendix B, we have encoded the discontinuity equations (3.1)- (3.4) in the simple statement that the P functions live on a two-sheeted Riemann surface and satisfy the relations (4.14), (4.15). However, some information is still missing to close the system. In fact, notice that, when describing the system with the functions P H ± , there is nothing to guarantee the correct large z behaviour (3.5) for the Y n,1 's with n ≥ 1. To fill this gap we will extract an extra equation from the TBA, equation (C.6), and show that it fixes the form of the gauge transformation (4.8)-(4.10). Finally, we will show that the final set of equations (4.14)-(4.16) is self-consistent and that equation (C.6) can be rederived from them.

Direct proof
The Y function sitting at the central node can be parametrised equivalently using the gauges T H or T V . This gives the two expressions

D. Dictionary
In this Appendix, we provide a dictionary between the notation of this paper and the ones of [16]. First, the conventions for the coupling constant and magnetic field are related as u here = U ref [16] /4, B here = H ref [16] /2 (D.1) The Quantum Transfer Matrix of [16] is characterised by three parameters: the Trotter numberN ∈ N, the inhomogeneity u and the rapidity v. The Gibbs free energy at temperature T can be obtained by taking the Trotter limit where Λ is the largest eigenvalue of the QTM. The map between v and the spectral parameter used in this paper, s, is 14 In particular, v = 0 corresponds to s → ∞.
Considering the identifications discussed in Section 5,