Semidefinite Programs at Finite Fermion Density

Semidefinite programs can be constructed to provide a non-perturbative view of the zero-temperature behavior of quantum systems. This paper examines the properties of these semidefinite programs when applied to lattice-regulated field theories exhibiting fermion sign problems. Specifically on the finite-density Thirring model, there is no indication that the accuracy of semidefinite programs suffers from any difficulty analogous to the sign problem.


I. INTRODUCTION
Monte Carlo methods, in the framework of lattice field theory, are often the only technique available for nonperturbative computations of expectation values in quantum field theories. These methods generally proceed by transforming a quantum expectation value to an equivalent classical expectation value in one higher space-time dimension. Unfortunately, for a wide variety of systems (notably including relativistic fermions at finite density and the Hubbard model away from half-filling), the resulting classical expression possesses a Boltzmann factor which is generally complex, and therefore cannot be viewed as a probability distribution for Monte Carlo sampling. Under standard complexity-theoretic assumptions, a theorem of Troyer and Wiese states that there is no fully general solution to these fermion sign problems [1].
Despite the nonexistence of a general solution to fermion sign problems, individual systems may be analyzed with a variety of approaches: complex Langevin [2], the density of states method [3], canonical methods [4,5], reweighting methods [6], series expansions in the chemical potential [7], fermion bags [8], analytic continuation from imaginary chemical potentials [9], and contour deformation methods [10]. These methods are frequently successful for low-dimensional (typically 1 + 1-dimensional) systems, but none have yet been able to provide reliable information about e.g. cold, dense quark matter.
A separate family of approaches to the zero-temperature behavior of quantum systems has recently been developed: the quantum-mechanical bootstrap [11][12][13][14]. Here we use the observation that, for any quantum mechanical operator O, the expectation value Tr ρO † O must be non-negative in any state ρ. As a direct consequence, if we select a basis of operators {O i } and a particular state ρ, we may construct a matrix guaranteed to be positive-semidefinite 1 : (1) * scott.lawrence-1@colorado.edu 1 Here and throughout, the notation M 0 is used to denote that a matrix is positive-semidefinite.
That is, for any vector v, we have v T M v ≥ 0. When the chosen basis of operators is complete-when any operator in the Hilbert space can be obtained as a linear combination of the elements of the basis-any set of expectation values consistent with Eq. (1) and the system's commutation relations is obtainable by some quantum state. In other words, every possible inequality describing the quantum system is captured by Eq. (1).
This bound by itself is not usually sufficient to provide interesting information about the ground state: after all, high-temperature states (and indeed far-from-equilibrium states) necessarily obey this bound as well. To do better we must include information about the Hamiltonian of interest. For a chosen Hamiltonian H, the ground state |Ω is by definition the state that minimizes Ω|H|Ω ≡ E 0 . Because the above bound captures every possible inequality, the relation Tr ρH ≥ E 0 can be derived from that bound.
The computational strategy of the quantum-mechanical bootstrap is to use Eq. (1) to derive lower bounds on the expectation value H , which must be obeyed by all states. These are bounds on the ground-state energy, and as more operators are used in the construction of M , the bounds converge to the true ground state energy. Concretely, we can cast the task of finding the ground state energy as an optimization problem over the space of expectation values, as follows. For a given basis of operators O i , we consider all expectation values of the form O † i O j . This defines the positive semidefinite matrix M ij of Eq. (1). The commutation relations furnish a set of linear constraints on the matrix elements of M -that is, a sequence of matrices C (k) such that we require We also impose that the expectation value of the identity operator is 1. This is equivalent to enforcing the correct normalization of the quantum states being looked at. Note that all of these constraints are linear in the expectation values, which is the space in which our optimization problem takes place. Finally, as long as the Hamiltonian is obtainable by a linear combination of our chosen basis, there is a matrix G such that Tr GM represents the expectation value of arXiv:2211.08874v3 [hep-lat] 9 Jul 2023 the Hamiltonian. We seek to minimize this linear combination of elements of M subject to the above constraints. Because this optimization problem is convex, the minimum obtained is guaranteed to be the global minimum, and therefore a lower bound on H that is valid for all states.
This type of optimization problem is termed a semidefinite program. This method provides a rigorous lower bound on the ground state energy. As a result, this method is essentially dual to variational methods, which provide similarly rigorous upper bounds on the ground state energy. The combination is particularly powerful when applied to systems for which the ground state energy is a quantity of central physical concern-for example in the study of quantum-mechanical bound states.
Separate from the explicit bound on the ground-state energy, we also obtain estimates of all expectation values used in the construction of the SDP. In situations where the ground-state energy is not of direct physical interest (most relevantly here, lattice-regularized field theories), these estimates are central to the utility of SDP-based methods.
When the basis of operators is complete, performing this minimization requires manipulating objects of size exponential in the volume of the physical system being studied. The same is true of attempting to exactly diagonalize the Hamiltonian, and so SDP-based methods appear to have no advantage. However, in practice it has been found that operator bases that are far from complete nevertheless result in remarkably tight bounds on simple quantum systems. This was reported early on for fewmatrix mechanics [11] and various one-body systems [12]. The power of the quantum-mechanical bootstrap derives from this observation: it is often possible to construct a polynomial-sized truncated basis for which the bootstrapped bound and many inferred expectation values are close to their true values.
Semi-definite programs have played an important role in quantum field theory before the advent of the quantummechanical bootstrap being used in this work. Most prominently, the conformal bootstrap provides the most precise constraints available so far on the critical expo-nents of many conformal field theories [15][16][17] (see [18] for a review of these methods). Bootstrap methods applied to the S-matrix, historically important [19], have also been recently revived [20,21].
This work applies SDP-based methods to quantum field theories by first regularizing the field theories on a spatial lattice (in the Hamiltonian formalism). The resulting system is a (many-body) quantum-mechanical system to which the methods of [11][12][13] are directly applicable. A similar approach was first used in [22] to treat the Hubbard model, and demonstrated directly on field theories in [23]. The hierarchy of positive-definiteness constraints so central to all quantum mechanical bootstrap methods is known as the Navascues-Pironio-Acin hierarchy [24,25].
The remainder of the paper is organized as follows. Section II introduces the Thirring model, a frequent test bed of methods targetting the fermion sign problem. The algorithms used to construct semi-definite programs, and extract from them expectation values, are detailed in Section III. Because the SDPs related to a space of expectation values, they can be constructed either at finite volumes or directly at infinite volume. The convergence behavior at finite volume, where direct diagonalization of the Hamiltonian is possible for comparison purposes, is examined in Section IV. Results at infinite volume are presented in Section V. Finally, Section VI discusses various aspects of the performance of this method, lays out open questions, and suggests avenues for future work.
All code used for this paper is available online [26]. SDPs are solved using the MOSEK toolkit [27].

II. THIRRING MODEL
The Thirring model is one of many afflicted by the fermion sign problem. Variants of this model have long been used as testbeds for methods to tackle sign problems [2,[28][29][30][31][32][33]. This work is in keeping with that tradition.
A lattice Thirring model with staggered fermions is defined by the Hamiltonian Here the second sum is taken over all pairs of adjacent sites on an L-site lattice with periodic boundary conditions. The operators c(x) are defined to obey anticommutation relations according to The lattice parameters m, µ, and g 2 are the mass, chemical potential, and coupling respectively. In the continuum limit, this model has emergent Lorentz invariance, and describes a single two-component fermionic field.
For positive g 2 , this model has a repulsive interaction between fermions (attractive between a fermion and antifermion), producing a fermion-antifermion bound state analogous to mesons in quantum chromodynamics. The mass ratio m B m F of the boson to the meson may be taken as a measure of the strength of the coupling-at weak coupling m B ∼ 2m F , and studies of the strong-coupling sign problem typically aim for m B ∼ m F [34]. This ratio m B m F is shown in Figure 1 across a range of couplings, with the fermion mass always tuned to be m F = 0.2 in lattice units, and the calculation performed by exact diagonalization of the Hamiltonian on a 10-site lattice. Calculations in this paper will typically be performed at g 2 = 0.5, well into the strong-coupling regime.

III. ALGORITHM
The basic algorithm to study a quantum system via a semidefinite program proceeds in four steps. First, a basis of operators {O i } is chosen. We will study the quantum system only by looking at expectation values of the form O † i O j ; all other properties are to be ignored. The axioms of our quantum system imply two sets of facts about these expectation values. First, the positive-definiteness of the inner product on the Hilbert space implies that the matrix M ij ≡ O † i O j is positive semi-definite. Second, the commutations relations imply certain linear relations between the expectation values. For example, for any sites x, y, we have c † (x)c(y) = δ xy − c(y)c † (x) . In the second step of the algorithm, these two sets of constraints are combined to construct a semi-definite program. The quantity to be minimized is the expectation value of the Hamiltonian, which can be expressed as a linear combination of matrix elements of M . The constraints are the linear relations between elements of M derived from commutation relations, and the nonlinear requirement M 0.
Once the SDP has been constructed, the third step of the algorithm is to perform the constrained minimization. For this work, the minimization is done using the MOSEK toolkit [27], which implements an interior-point method to perform the optimization [35].
Finally, in the fourth step we wish to extract physical information from the solution of the SDP. For a simple quantum-mechanical system, the lower bound on the ground state is of intrinsic interest. In the context of (lattice-regularized) field theories, the ground state energy density at a single value of µ has no physically meaningful continuum limit. For zero-temperature, finite-chemicalpotential models, the observable of greatest interest is typically the number density n(µ). Alternatively, the difference between ground-state energies 2 at two different values of µ provides the integral of n(µ) according to When solving an SDP, typical algorithms (including the interior-point methods used by MOSEK) return not only the optimal value of the objective function (in this case corresponding to the ground-state energy), but also a value of the positive semi-definite matrix M that obtains that optimum. Any observable encoded in the matrix M can therefore be directly estimated once the optimization is complete. In the case of the Thirring model, because a term µN appears in the Hamiltonian, the fermion number will always be accessible from the expectation values present in M .
The fermion number could also be obtained, in principle, from solving two SDPs for slightly different values of the chemical potential µ, and computing N = − ∂ ∂µ Ĥ − µN . In fact this yields the same result as simply reading off the value of N from the matrix M in the solution to the SDP. To see this, consider the positive semi-definite matrix M obtained after solution. Generically, this matrix will have exactly one vanishing eigenvalue. The corresponding eigenvector v defines some operatorB ≡ v i O i such that the expectation value B †B is estimated to be 0. Using the linear part of the SDP (i.e. the commutation relations), it may be shown that B †B =Ĥ − E, where E is the SDP-obtained estimate of the ground-state energy. Now let us consider what happens after perturbing the Hamiltonian by a term N . Again generically,B is also perturbed by some term of the same orderB →B + Ĉ . We now see that (expanding to first order in ) where again the equalityB †Ĉ +Ĉ †B =N follows from the commutation relations. As a result, the change in the estimated ground state energy is equal to the expectation value of this operator, which is constrained by the commutation relations to be equal to the estimate N . Typical computational methods (e.g. lattice Monte Carlo methods) arrive at the infinite-volume limit by performing calculations at a sequence of finite volumes and then extrapolating to obtain the infinite-volume result. As discussed in [23], because SDP-based methods work directly with operators and make no reference to the Hilbert space itself, it is natural to perform calculations directly in the infinite-volume limit. This is accomplished with two modifications to the algorithm described above. First, instead of minimizing the Hamiltonian expectation value (which is no longer well-defined), we minimize the Hamiltonian density. Second, we impose translational invariance 3 by including in the formulation of the SDP additional linear constraints of the form O = T † OT for any unitary translation operator T .
The core of the quantum-mechanical bootstrap is the choice of an incomplete basis of operators-in many ways this parallels the need for a choice of basis when explicitly diagonalizing a truncated Hamiltonian. The time required to solve an SDP scales roughly as the cube of the size of the basis, so it is critical to keep the basis small. At the same time, a poor choice of basis (including an insufficiently large basis) will not yield a tight bound on the ground-state energy. Finally note that, although the approximation to E 0 obtained via a hierarchy of SDPs is systematically improveable in the sense that, with a sufficiently large basis, the approximated E 0 can be brought within any desired tolerance of the exact result, in principle the cost may be exponential in the precision requested. An efficient algorithm depends on a well-chosen basis.
For this paper, we define a hierarchy of four bases. Each basis is both translationally invariant and Hermitian (meaning that for any operator O in the basis, O † is also present). The first, labelled H0, consists of c(x), c † (x), and c † (x)c(x) for all lattice sites x. On an L-site lattice, counting the identity, this has 3L + 1 elements. Note that although the Hamiltonian is not in the span of this basis, it is in the span of operators obtained by taking a product O † i O j of two elements of this basis, and so the SDP is able to provide a nontrivial constraint on the ground-state energy.
The next basis, H1, adds to H0 the L operators of the form c(x)c † (x + 1), and their L complex conjugates. This basis has a total of 5L + 1 elements.
To the elements of H1 we can add the L operators found in the interaction term of the Hamiltonian: for each site x, c † (x)c(x)c(x + 1)c † (x + 1). The basis obtained in this way is labelled H2.
These first three bases are of course inspired by terms found in the Hamiltonian of Eq. (3). The final basis, termed C1, adds to H2 the L operators of the form c † (x)c(x + 1)c † (x + 2)c(x + 2), and their complex conjugates. The size of this largest basis is 8L + 1.
To check the correctness of these methods, Figure 2 provides a comparison on a 10-site lattice with ground state energy information obtained by exact diagonalization of the Hamiltonian. The left panel shows the ground-state energy bound proven via SDP. Although the bound itself is not meaningful, an estimate of the integral of number density between two chemical potentials may be obtained from the difference of energies at those values of µ. (The extra term proportional to µ ensures that the Fermi sea is not counted in the estimated fermion number.) Alternatively, the right panel compares the estimated fermion number to the exact result. For all four bases, the energy bounds are seen to lie below (or on) the exact value, as they should.
Moreover, significant convergence towards the true value is seen as the basis is expanded. On this lattice, the largest basis (C1) has 81 elements. By comparison, a complete basis of operators on the 1024 states would require 2 20 operators.
The next section will consider the convergence properties of the SDPs obtained from these bases more carefully. Already from Figure 2, we can see that the basis H0 provides a qualitatively accurate equation of state everywhere away from the step at µ ∼ 1.8. The mechanism for this is visible in the left panel of that figure: for µ ∈ [0, 1.5), the error in the energy is roughly constant.

IV. CONVERGENCE
It has previously been demonstrated that semidefinite programs can yield efficient and accurate approximations for the ground-state properties of quantum lattices [22,23]. Other methods, notably lattice Monte Carlo methods, suffer at finite chemical potentials and strong couplings due to the fermion sign problem. The purpose of this section is to investigate to what extent the same is true of the SDPs defined in the previous section. Figure 3 shows the error in the ground-state energy estimate for all four bases described above. The most striking feature of this figure is that all SDPs perform best at large chemical potentials-this is the opposite of the 'traditional' behavior of Monte Carlo methods. This success is due to the fact that at large chemical potential, the true ground state of the system is saturated, with all fermionic modes being occupied. Any SDP will, at sufficiently large µ, prefer this state to all others and therefore agree with the exact result.
As more operators are added to the basis, the SDP approximation improves. Again, the largest and easiest improvements are achieved at larger chemical potentials. In fact, of the bases considered, only the largest (C1) is substantially different at µ = 0.
The news is not all good. Early attempts to resolve the fermion sign problem of lattice QCD and related models frequently struggled to replicate the so-called silver blaze phenomenon [36]. (Contour deformations methods-such as those based on Lefschetz thimbles-were a notable exception, e.g. [37].) In the infinite-volume limit, this term refers to the fact that the number density remains constant for some interval of chemical potentials µ ∈ [0, µ c ) before suddenly beginning to increase. In QCD, this increase begins slightly below the proton mass-faulty algorithms for alleviating the sign problem would characteristically begin to increase far too early. We see that the SDPs investigated here mostly (H2 is a striking exception) have a similar flaw: the density begins to increase near the bare fermion mass of 0.05, rather than the correct physical mass m F ≈ 0.20.
Moreover, at a finite volume, all increases in the fermion number ought to be discontinuous jumps: the fermion number in a non-degenerate eigenstate of the Hamiltonian is always an integer. The semidefinite programs do not generally have this behavior, instead appearing to smooth out these jumps. Again, the basis H2 is a striking exception, indicating that there may be some condition on the basis which is sufficient to enforce this desired behavior.
The performance of the SDPs as coupling is changed is more consistent. The left panel of Figure 3 shows the error in the estimate of the ground-state energy at µ = 0.5, for a sequence of theories of increasing coupling g 2 . All computations are done on a 10-site lattice, with the bare mass m tuned such that the renormalized fermion mass (measured on the same lattice) is m F = 0.2. As the coupling is increased, the SDP's error increases roughly linearly.
All bases considered give an exact bound for the free theory. For this theory, the raising and lowering operators that diagonalize the Hamiltonian can be constructed via the discrete Fourier transform of the position-space creation and annihilation operators c † (x), c(x). Therefore, these operators diagonalizing the Hamiltonian are in the span even of the smallest basis, H0. This is sufficient for the SDP estimate to be exact.

V. INFINITE VOLUME
As noted in [22], because the semidefinite programs constraining expectation values work in terms of opera-tors rather than configurations or states in Hilbert space, these calculations can be performed directly in the infinite volume limit. Of course, even for a system (like the Thirring model) with a finite-dimensional local Hilbert space, the set of operators available in this limit is infinite. As a result, any finite basis chosen to construct the SDP will yield only an approximation. Nevertheless, each estimate of the ground-state energy density will be a rigorous lower bound to the infinite-volume result, and we may hope that the estimates of the number density provide a good approximation in practice.
In the infinite-volume limit, it is no longer possible to speak of the expectation value of the Hamiltonian. Instead, we must work with a Hamiltonian densityĥ(0), defined so that the Hamiltonian isĤ = rĥ (2r): The optimization objective for the SDP will be the expectation value ĥ (0) . In order to close out the SDP, we must include some information about how the energy density around site x = 0 is related to the energy density elsewhere on the lattice. In this case we will assume translational invariance. In particular, for any operator O and even displacement 2r (r ∈ Z), we assume that All linear constraints of this form are added to the SDP. At finite volume, the free theory is exactly solved by an SDP constructed from a basis consisting only of all available creation and annihilation operators. At infinite volume that basis already must be truncated in order to obtain a finite SDP. Any such truncation will have some maximum radius of operators L, and will be an approximation that improves as L is taken to be larger.
The difference between the truncated SDP and the exact result can be thought of as a finite-volume effect. Consider an SDP constructed from the operators a(−L), . . . , a(L) and their complex conjugates-or, in the case of a theory at non-zero coupling, an SDP consisting of all possible operators with support only on [−L, L]. Such an SDP captures all quantum mechanical constraints on expectation values that arise on that range, but has no information about the boundary conditions. As a result, the obtained bound on the energy density can only be as tight as the lowest energy density achievable with any boundary conditions.
The performance of ever-larger SDPs for the free theory is shown in Figure 4. Each SDP is constructed only from the operators a(−L), . . . , a(L + 1), and their complex conjugates. The SDP-estimated number density is compared to the exact lattice result The most notable feature is that the approximation converges much faster, as a function of L, when the fermion mass is larger. Physically, this is a result of the fact that finite volume effects decay parameterically with the dimensionless mL. The basis used to construct the SDP, described in the text, grows linearly with the parameter L. For all bare masses the SDP converges to the exact number density, but the convergence is somewhat faster for larger masses, where only short-distance correlations are relevant. All calculations are performed at µ m = 4 3 . To make the differences in convergence more apparent, the bottom panel shows the relative error on a log scale.
To perform calculations with a nonzero interaction, it is necessary to return to larger bases of the sort considered in the previous two sections. Infinite-volume equivalents of the bases H0, H1, H2, and C1 are defined by simply requiring −L ≤ x ≤ L + 1. Figure 5 shows calculations performed in the infinite volume limit for interacting theories. The left panel shows a weakly coupled theory (g 2 = 0.1), and the right a strongly coupled theory (g 2 = 0.5).
In both theories, the number of visible plateaus increases as the permitted separation L of operators is increased. This again shows the manner in which the finite-L effects are analogous to finite-volume effects: calculations at finite volume will also show a number of plateaus proportional to the volume of the system.
The weakly-coupled theory in particular shows good qualitative agreement between the various approximations. Although no extrapolation is performed here, both show signs of convergence as larger bases are considered.

VI. DISCUSSION
This paper has explored the performance of SDP-based methods in analyzing the Thirring model in two spacetime dimensions, at finite fermion density. The SDPs used are small relative to the size of the Hilbert space of the system under study, and there is no indication that the method encounters any obstacles analogous to a fermion sign problem. While Monte Carlo methods perform best at µ = 0 and worst when µ is large enough for the lattice to be saturated, the SDPs have the reverse behavior, with tight bounds most easily achieved for chemical potentials at or near lattice saturation density.
Fermion sign problems are most often alleviated while remaining within a Monte Carlo framework. It is difficult to construct a fair performance comparison between SDPs and these Monte Carlo-based methods. For example, if we consider the asymptotics of the algorithm when the answer is required to be accurate to within , we see that SDP methods formally outperform Monte Carlo methods. For a system with a finite-dimensional Hilbert space, a Monte Carlo method requires time exponential in the requested number of significant figures. Since, for a fixed finite-dimensional system, there is an SDP that obtains the exact answer, the time required for the SDP to perform a computation to within that precision is of order 1, formally 4 -but that observation provides us with no practical information about the performance of SDPs.
The SDPs considered here provide qualitatively accurate equations of state at finite volumes, across the entire range of chemical potentials. However, the equations of state obtained fail a few important tests. Most notably, the estimated number density, and the "silver blaze" phenomenon at small chemical potentials is not correctly reproduced.
As with other methods along the lines of the "quantummechanical bootstrap" [11][12][13][14]23], the method in this paper is in a sense dual to the usual variational principle. Whereas a variational study of a quantum-mechanical system yields a rigorous upper bound on the energy, the solution to an SDP yields a rigorous lower bound to the