The SU(3) spin model with chemical potential by series expansion techniques

The $SU(3)$ spin model with chemical potential corresponds to a simplified version of QCD with static quarks in the strong coupling regime. It has been studied previously as a testing ground for new methods aiming to overcome the sign problem of lattice QCD. In this work we show that the equation of state and the phase structure of the model can be determined to reasonable accuracy by a linked cluster expansion. In particular, we compute the free energy to 14-th order in the nearest neighbour coupling. The resulting predictions for the equation of state and the location of the critical end point agree with numerical determinations to ${\cal O}(1\%)$ and ${\cal O}(10\%)$, respectively. While the accuracy for the critical couplings is still limited at the current series depth, the approach is equally applicable at zero and non-zero chemical potential, as well as to effective QCD Hamiltonians obtained by strong coupling and hopping expansions.


Introduction
Despite growing demand from various fields of physics, the details of the QCD phase diagram as a function of temperature T and baryon chemical potential µ B remain unknown to date. This is because of a severe sign problem due to the complex fermion determinant for µ B = 0, which prohibits straightforward Monte Carlo simulations of lattice QCD (for introductions, see, e.g., [1,2]). Controlled results only exist by indirect methods for the low density region µ B < ∼ 3T , where no sign of criticality is observed [3,4].
This has motivated the search for alternative formulations and algorithms to solve this problem. While there is a vast literature on the general subject of sign problems, all approaches devised so far work for limited classes of Hamiltonians, which do not (yet) include QCD. The purpose of the present paper is to test series expansion techniques, generically known as "high temperature" expansions in the condensed matter literature, against the known numerical results for the SU (3) spin model.
The SU (3) spin model with chemical potential is a system often used in the literature to test new methods aiming at the QCD sign problem. It corresponds to QCD with a simplified determinant for static quarks near the strong coupling limit. It features a local SU (3) symmetry whose center Z(3) spontaneously breaks at some critical coupling, as well as a sign problem at µ = 0. The model was reformulated free of a sign problem in terms of a flux representation [5], allowing for simulations [6] by means of a worm algorithm.
Likewise, simulations using a complex Langevin algorithm have been successful [7][8][9]. The phase structure of the SU (3) spin model is therefore known at zero and finite chemical potential with good precision.
Being analytic in nature, series expansion techniques are insensitive to sign problems and in principle applicable to any Hamiltonian resembling a spin model. This is of particular relevance to finite density QCD, where combined strong coupling and hopping expansions produce effective Polyakov loop actions which closely resemble the model considered here [10][11][12][13]. If such series expansions can be pushed to sufficiently high orders for a required accuracy, the extension to µ B = 0 poses no fundamentally new problem. In this work we demonstrate that for the SU (3) spin model the equation of state as well as the phase diagram can be determined to satisfactory accuracy with a reasonable effort.
The paper is organised as follows. Section 2 introduces the spin model as well as the observables which we use to characterise its phase diagram. In section 3 we briefly summarise the main features of the linked cluster expansion and apply it to the spin model at hand. From the resulting power series, we then extract the equation of state and compare with Monte Carlo results at zero density. We discuss how to extract the phase structure and its critical features by means of series analyses in section 4. Finally, section 5 concludes with a discussion of the prospects of this approach for effective theories of QCD.

The SU (3) spin model
We consider the action τ L(x)L * (x +k) + L * (x)L(x +k) + κ e µ L(x) + e −µ L * (x) , (2.1) with x denoting the sites of a three-dimensional cubic lattice withk the unit vectors in the three directions. The fields L(x) = TrW (x) are complex scalar variables, resulting by taking the trace of the SU (3)-matrices W (x). The analogy to other spin models is emphasised by defining Then L plays the role of magnetisation and η,η represent symmetry breaking external fields. For η =η = 0 the action is invariant under global transformations It is obvious that for µ = 0 the action is real. The partition function then is straightforward to simulate by standard Metropolis methods, which we will use as a benchmark for our series expansions. In this case we work with finite volumes and periodic boundary conditions in all directions. For µ = 0, the action is complex and the simulations exhibit a sign problem.
The natural observable to compute by series expansion is the free energy density which is related to the trace of the energy momentum tensor and in QCD is called interaction measure. In order to identify phase transitions, we look for maximal fluctuations in the spin model analogues of the magnetic susceptibility and the specific heat, respectively,

Linked cluster expansion
There are many expansions in physics that can be expressed in terms of linked clusters (connected graphs). For an introduction to the general subject see [14], where the method we are employing is also called the 'free graph expansion', a name which will become clear later. One famous application of this expansion is to the φ 4 theory [15], where it was used to prove triviality [16]. Extensions to very high orders in this case can be found in [17]. Following these references, we employ a so-called 'vertex-renormalisation' 1 , where the number of graphs one has to consider is decreased at the cost of increased algebraic complexity. Further improvements can be achieved by an additional 'bond-renormalisation'. An impressive example of applying both renormalisation techniques is [18], where two-point Green's functions for a generalized 3D Ising model are obtained to 25th order on the bcc lattice and to 23rd order on the sc lattice. Our application follows closely the review article by Wortis [19], which gives a clear explanation of the method.

Graphs
In this section we review what is necessary to give a precise meaning to the involved formulas. Our notation and definitions are based on a combination of those introduced in [16,18,20,21]. One complication we encounter is that, due to the particular form of the interaction term in the spin model, we have to consider directed graphs. A directed graph consists of two finite sets, a set V (G) of vertices and a set D(G) of directed bonds together with two mappings In this way, each bond b ∈ D(G) is associated with the ordered pair (i(b), t(b)) of endpoints, where i(b) is called the initial and t(b) is called the terminal point of b. The bond can be visualised as an arrow with tail i(b) and head t(b). Two vertices, which are the endpoints of the same bond, are adjacent to each other. The in-degree n in (v) of a vertex is the number of directed bonds with t(b) = v and the out-degree n out (v) is defined analogously. Their sum gives the degree of a vertex, Sometimes one distinguishes external/rooted vertices R(G) and internal vertices I(G), V (G) = R(G) ∪ I(G). An r-rooted graph denotes a graph with r external vertices and we will always assume that the external vertices have the labels 1, . . . , r.
On the level of graphs, the external and internal vertices differ in how they enter the definition of a graph isomorphism. Two graphs G and G are considered isomorphic if there exist two mappings with the properties that (a) ϕ and λ are bijective, From now on, whenever we define a set S of graphs, we always actually mean the quotient set S/ ∼, where G ∼ G if they are isomorphic. The symmetry factor S(G) of a graph G is the order of the automorphism group of G. Two vertices v 1 and v 2 are called connected, if there is a sequence of adjacent vertices w 1 , . . . , w n , with w 1 = v 1 and w n = v 2 . A graph is defined to be connected if its vertices are pairwise connected.
We write G \ v for the deletion of the vertex v from the graph G. This means that For a connected graph with external vertices, an articulation point is a vertex v with the property that G \ v has a vertex which is not connected to a root. In the case of a graph with no external vertices, an articulation point is a vertex v such that G \ v is disconnected. A 1-irreducible graph is then a graph with no articulation points.
Finally, a 1-insertion denotes a connected graph with one external vertex, where deleting the external vertex leaves the graph connected. In case of a 1-insertion we define 1-irreducibility to mean that there is no articulation point except of the root (otherwise, the 1-irreducible 1-insertions would only consist of one single vertex graph).

Unrenormalised expansion
We use the linked cluster expansion to determine the Taylor expansion of the free energy of the spin model in τ around τ = 0 in terms of connected graphs. The general graphs which enter this expansion are valid for many different theories, so one has to make a connection between the graphs and the concrete theory. This is achieved by assigning weights W 0 (G) to the graphs. These weights are expressed in so-called bare semi-invariants, which are the cumulant correlations of the theory without interactions.
In the case of the SU (3) spin model the bare semi-invariants M 0 n,m (η,η) with n, m ∈ N, It is useful to have an expansion of the logarithm in terms of η andη. To this end, consider Using Faà di Bruno's formula one then obtains where B n,k denotes the partial Bell polynomial (3.8) We use the notation that any logical statement in square brackets evaluates to 1 if it is true and 0 otherwise. The derivatives of z can be computed via for which an explicit formula is given in [5]. For a directed graph G we can now define its weight W 0 (G) to be So far, all definitions are independent of the underlying lattice of the theory. In the linked cluster expansion, the information about the lattice is encoded in the free embedding numbers, also called free multiplicities. On a three-dimensional cubic lattice, the free embedding number counts the number of ways to put a graph on Z 3 , with one vertex fixed to the origin while for all others adjacent vertices have to correspond to nearest neighbours on the lattice. For the free embedding number it is allowed to place two vertices on the same lattice point, which makes its computation relatively easy in comparison to other embedding numbers. To give a more formal definition, let X G,v denote the set of all functions that map all vertices of the graph G to where d denotes the lattice distance We are now able to write the expansion of the free energy density in terms of connected graphs. Define G (i) c to be the set of all connected and directed graphs with i bonds and no external vertices, then the free energy density evaluates to To be a bit more explicit, we evaluate this formula to O(τ 2 ) as an example. The involved graphs, their symmetry and embedding numbers and their weights are listed in table 1.
Consequently, one has The bare semi-invariants can then be evaluated according to the explanation at the beginning of this section.

Vertex-renormalised expansion
The set G (i) c becomes large rather fast when i is increased. As mentioned above, here we employ a method called vertex-renormalisation, which reduces the number of graphs one has to consider.
To this end, we first define the self fields, which collect all possible ways to decorate a vertex (note that the free embedding number factorises along articulation points). Denoting by I (i) (n,m) all one-insertions with i edges and where the external vertex has out-degree n and in-degree m, the self field G n,m is defined to be (remember that W 0 was defined in a way such that external vertices give no contribution) This in turn can be used to establish the notion of renormalised semi-invariants (3.17) Defining the renormalised weight W (G) of a graph in the obvious way, one can obtain the free energy density via contains sums over G (i) 1 , the sets of all 1-irreducible graphs with no external vertices and i edges. Those sets are much smaller than G (i) c . At this point, however, nothing is gained yet, because in order to obtain the self fields one still has to evaluate sums over all 1-insertions (also those which are not 1-irreducible) in equation (3.16). Replacing all bare semi-invariants by their renormalised counter parts in this equation enables one to restrict the sums to 1-irreducible 1-insertions I As a result, however, equations (3.17) and (3.20) are coupled equations. Therefore, the following iterative strategy is used: suppose one wants to obtain the free energy to the order ∼ τ nmax . Using equation (3.18) this necessitates the determination of the renormalised semi-invariants M n,m to order n max − n − m. To leading order which enables the determination of G 1,0 and G 0,1 to order 1 in τ using equation (3.20). These self fields can then be used to determine the renormalised semi-invariants to order 1, and in this way the procedure can be continued to the necessary order. In general, having obtained all renormalised semi-invariants M n,m to order p with (n + m) ≤ p, one can determine the G n,m to order p + 1 with (n + m) ≤ p + 1 and use those to obtain the renormalised semi-invariants to order p + 1.
We again illustrate the procedure by deriving the free energy to order O(τ 2 ). At first, we have to determine the renormalised semi-invariants M n,m and the self-fields G n,m to O(τ 2 ). Since the leading contribution of G n,m is of O(τ n+m ) the relevant contributions from equation (3.17) to the renormalised semi-invariants are Furthermore, the relevant 1-irreducible 1-insertions for the self-fields are shown in table 2 and therefore  This in turn can be used in equation (3.22) to obtain Inserting these equations into equations (3.23) to (3.25) results in Next one determines the Φ-functional from equation (3.19). The relevant graphs are graphs 2 to 4 in table 1, however with the bare semi-invariants replaced by their renormalised counterparts for the weights: Inserting the expressions for the renormalised semi-invariants and the self-fields in terms of the bare semi-invariants which were derived above into equation (3.18) then results in the same expression that was obtained in equation (3.15). While at this point it might seem that the renormalised procedure introduces unnecessary complications, we remark that for higher orders a significant reduction of graphs is achieved.
redundancies. To this end, we introduce the following order-relation on multi-indices of N k : for n, m ∈ N k one has where the combinatorial factor m(l 1 , . . . , l p ) counts the number of unique permutations of its arguments. For graph isomorphism checking and symmetry numbers we used the graph procedures provided by McKay's nauty [22]. The embedding numbers were computed in Mathematica using the algorithm described in [16]. The computation of the graph weights and the vertexrenormalisation were also done in Mathematica, among other things because it implements multi precision arithmetic via the GMP [23] in a convenient way.

The equation of state
Using these tools, we derived the free energy through O(τ 14 ) and, in order to keep the size of the expressions under control, all products of bare semi-invariants were expanded and The equation of state is often given in terms of the so-called interaction measure equation (2.6), since the latter can be straightforwardly determined in Monte Carlo simulations. All thermodynamic functions can be computed from it by integration or differentiation. Our series results to the three highest orders are shown in figure 1, where good convergence is observed. At vanishing chemical potential (left), we compare our series directly to Monte Carlo data. We find excellent quantitative agreement up to the immediate neighbourhood of the phase transition, which a finite polynomial can only indicate by loss of convergence. On the right we show results for a large chemical potential, which poses no fundamental problem for a series expansion method, whereas standard Monte Carlo simulations are no longer possible because of the sign problem.

Resummation by Padé approximants
Finite series generally break down in the vicinity of phase transitions. A marked improvement in convergence properties can often be obtained by infinite-order resummations, like mappings of the expansion variables, use of renormalisation group techniques or approximation by rational functions. For a general review and introduction, see [24]. Here we model a function f (x), known only as finite power series, by Padé approximants defined as rational functions,

The phase transition
In the previous section we have seen that the linked cluster expansion gives excellent results for the equation of state in the symmetric phase of the SU (3) spin model at zero as well as finite chemical potential. In this section we explore various methods to extract information about the location and order of the phase transition. Our preferred observable to locate phase transitions is the susceptibility, equation Because of its definition through τ -derivatives, the series for the specific heat, equation (2.8), is by two orders shorter and, moreover, not related to the order parameter. We observe it to generally exhibit poorer convergence near the transition.

Radius of convergence from the ratio test and Padé approximants
When a function with a given domain of analyticity in its complex argument is expanded in a power series, the radius of convergence is given by the distance between the expansion point and the nearest singularity. If this singularity is on the real positive axis, it signals a phase transition. There are various estimators for the radius of convergence, which are appropriate for different circumstances. The simplest and most well-known is the ratio test of consecutive coefficients of a series, which is expected to converge whenever coefficients have either the same or alternating signs. The radius of convergence is then obtained by extrapolation, r = lim n→∞ r n , r n = c n c n+1 .
which suggests a linear extrapolation to determine the radius of convergence. (For sufficiently long series and second-order transitions, λ is related to a critical exponent. Here we merely treat it as a fit parameter.) We indeed observe a linear trend for the higher ratios for many points in the parameter space, but with a superimposed oscillatory behaviour. This is familiar from Ising models on loose-packed lattices (in that case due to an antiferromagnetic singularity) and can be reduced by using the shifted ratios [25] (4.4) Figure 2 shows the shifted estimators r −1 n plotted vs. n −1 for two values of κ and µ = 0. As is apparent in the figure, a rather accurate estimate for the radius of convergence is obtained on the left, while on the right the linear scaling region has not yet been reached.
Repeating this procedure for varying values of κ, we obtain a line of singularities which is shown in figure 3 (left), together with the phase structure obtained by Monte Carlo simulations. The latter show a line of first order transitions, which weaken with κ to end in a critical point, which has been computed in simulations of a flux representation [6]. We observe the radius of convergence estimate from our series to systematically overshoot the critical coupling by an amount which diminishes towards the critical endpoint, where it vanishes. This behaviour is due to the fact that the series has information from one phase only. Thus the singularity detected by a series analysis is the end of the metastability region of that phase, rather than the true critical coupling, which requires information from both phases. The same behaviour is observed in, e.g., Potts models [26] with firstorder transitions.
A series analysis well-suited for second-order transitions is provided by Padé approximants. At a second-order transition, the susceptibility diverges with a critical exponent and, approaching the transition, its logarithmic derivative has a simple pole with the critical exponent as its residue, Such a pole in the Dlog can be faithfully reproduced by Padé approximants. A pole in the full function D χ is thus predicted whenever different approximants show converging roots of their denominators. While Padé approximants work astonishingly well for one-variable problems with a second-order transition (as a recent example, see the strong coupling expansion of SU (2) Yang-Mills theory [27,28]), the two-variable situation of the SU (3) spin model with different types of transitions is more complicated. Figure 3 (right) shows the results of these estimates, which are in remarkable agreement with those of the ratio test. The critical point is accurately reproduced at the appropriate κ-value, and the approximants also appear to pick up the end of the metastability range in the first-order region. However, both the ratio test and the Padé approximants indicate singularities also in the crossover region, where the full χ is known to be analytic. (With sufficiently long series, the radius of convergence should either become infinite or move to a complex value in this regime.) In summary, both analyses pick up the rise of the susceptibility near the phase boundary, but are unable to clearly distinguish between orders of the phase transition or crossover behaviour. Unfortunately, this difficulty remains when multi-variable generalisations of Padé approximants, so-called partial differential approximants (see e.g. [24]), are used in both variables, as we have explicitly tested.

The critical point
In order to locate the critical point, we now exploit the following facts: first, Dlog Padés model the singular behaviour of the full observables correctly only at a second order point; second, at such a critical point the susceptibility χ and the specific heat C show the same   Table 4. Intersection points of the poles in the Padé approximants to D C , D χ , together with the residues at the intersection.
with the same mixed exponent [29,30], For the 3d-Ising universality class, we have γ = 1.24, β = 0.326, δ = 4.79 and = 0.79. Figure 4 shows the location of the poles of the two highest-order approximants to each D χ and D C . According to the arguments given above, these need to agree at a second-order transition and only there, so we take their intersections as estimates for the location of the critical point. Note the comparatively worse convergence of the poles in C, which causes the largest part of the uncertainty.
We now systematise this analysis to obtain an error estimate. In table 4 we list all intersections of diagonal Padé approximants for D C and D χ as estimates for the critical point. The last two columns give the residue associated with each observable at the intersection point. Padé approximants necessarily produce ever more poles, the higher their order. It is then clear that some of them are artefacts and do not have physical meaning. To get rid of these, we demand that the residues of the intersecting observables agree within 20%. These estimates are then averaged over.
The upper part of table 4 corresponds to µ = 0 as in figures 3, 4. Averaging over the different estimates gives our final result for the location of the critical point, which is shown in figure 5 (left), in comparison to the simulation result from [6]. The same kind of analysis can be repeated for finite chemical potential without further difficulties.  [6]. The predictions from the series approach have a relative error of about 10% in τ c and about 20% in κ c or µ c . Presumably the larger error on the latter variables is due to the flatness of the critical line in the phase diagram, figure 3, so it takes more accuracy (and thus higher orders) to resolve changes in those directions. Within error bars, the predicted critical points agree with the numerically determined ones, so the estimate of the systematic error is realistic. Finally, we compute the critical exponent from the more reliable D χ and find = 0.69(2) and = 0.67(8) for the two computed critical points, respectively. This is within ∼15% of the true value.

Conclusions
We have studied the SU (3) spin model with chemical potential employing analytic series expansion techniques. In particular, we determined the free energy density through order τ 14 in the coupling of the energy-like term and through order κ 60 in the magnetisation-like coupling. In the phase with unbroken center symmetry, where the expansion is valid, we observe excellent convergence and the numerically known equation of state is reproduced to ∼1% up to the phase transition. The series holds for both zero and finite chemical potential, since it is unaffected by sign problems. We then investigated different approaches to extract information about the phase transition from the power series. Agreeing estimates for the radius of convergence were obtained by extrapolations of the ratio test and Padé approximants to the series of the susceptibility in the τ -variable. In the parameter region with a first-order transition, the detected singularity corresponds to the end of the metastability region of the symmetric phase. At the critical end point of the transition, the estimates agree with the true transition, but finite radius of convergence estimates continue far into the crossover region, where there is no singularity for real couplings. Such methods by themselves are therefore unable to determine the nature of a phase transition in a multi-variable system. On the other hand, exploiting the fact that the magnetic susceptibility and the specific heat diverge in the same manner at the critical point, an analysis of the crossings of their respective Dlog Padés allows to extract estimates for the location of the critical point with a relative accuracy of 10-20% in the critical couplings.
In comparison to numerical solutions via flux representations, the current depth series is competitive for the equation of state, but clearly less accurate for the critical couplings. Nevertheless, it offers an interesting perspective to finite density QCD, where such simulations are not yet applicable. In particular, our results show that effective theory approaches based on strong coupling and hopping expansions [10][11][12][13] are reliable and, with some effort to develop efficient computational schemes for higher orders and multiple couplings, can be systematically improved to apply to larger regions of the QCD parameter space.