Extension of the Uhlenbeck-Ford Model with an Attraction

The Uhlenbeck-Ford model for soft repulsion, which has only a repulsive interaction, is extended by inclusion of an attraction. This extension still allows an analytical evaluation of the virial coefficients. The integrals over the graph contributions are reduced to a combinatorial problem. We have calculated the virial coefficients to order 6 in the density. A link is made between this model and more common interactions, like the 12-6 Lennard-Jones potential.


Introduction
The expansion of the thermodynamic properties in terms of a power series in the density is the oldest tool in the systematic study of the properties of dense gases. It has the form [1] p where p is the pressure, T the absolute temperature and n the number density. The coefficients B l are called the virial coefficients (B 1 = 1). The expansion is unique in the sense that each term in the series has an explicit prescription for its calculation. Mayer [1,2] was the first to introduce a graphical representation for the various contributions. For low densities a few terms suffice and in this domain the virial expansion has proven to be an invaluable tool for the computation of the thermodynamic properties. Not only the pressure but also the other thermodynamic quantities can be expanded in a series in the density n expressable in the same coefficients B l . A similar expansion for the transport properties has been attempted [3,4,5,6], but even the first correction to ideal gas behavior leads to divergencies and it turns out that the transport properties are not expandable in a power series in the density. The virial coefficients B l are found as l-fold integrals B l = − l − 1 V l! dr 1 dr 2 · · · dr l H l (r 1 , · · · , r l ), where the function H l is represented by graphs with l vertices and a number of occupied edges, each carrying the (Mayer) function with V (r ij ) the intermolecular potential between the particles i and j. The problem of evaluating the virial coefficients is twofold. The first part is the generation of the graphs and the second part are the integrals associated with the graphs. The generation of the graphs is general and independent of the interaction potential of the gas, while the evaluation of the integrals is strongly dependent on the interaction. Usually finding the graphs is not the limiting factor since the integrals become already for order 5 a too-complicated for e.g. a Lennard-Jones interaction. An exception forms the so-called hard-sphere interaction where the virial coefficients are evaluated up to order 10. The hard-sphere model has been thoroughly investigated in all dimensions [7,8,9,10], in particular in its relation to the convergence of the virial series.
The generation of graphs has been extensively studied by the mathematicians, who have developed an efficient algorithm for this problem. Even with their efficiency the enumeration gets laborious around order 10, due to the more than exponential growth of the number of relevant graphs. For our purpose the bare generation of the graphs is not sufficient, we have to know also the symmetry properties of the graph. As this is a job valid for all systems, one has to generate the sequence only once and the time spend on it is generally not the bottle-neck for the calculation. By straightforwardly generating all graphs and eliminating those that do not qualify, we could reach all the graphs up to l = 9 vertices together with their symmetry properties.
It is an old idea, due to Uhlenbeck and Ford [11], to simplify the evaluation of the graphs by replacing the Mayer function f (r) by a gaussian One can view this replacement as an approximation of a potential, but also as a model on its own. As the corresponding potential is positive everywhere, with an infinite limit for r = 0, the Uhlenbeck-Ford (UF) model can be considered as representing a soft repulsive system with range a. The advantage of this model is that the integration of any graph contribution becomes a gaussian integral, which can be evaluated analytically. The model Eq. (4) contains only one parameter, the length scale a which can be combined with the density n to a dimensionless measure for the density.
Recently this idea has been vigorously picked up by Leite et al. [12,13]. They evaluated not only the virial coefficients up to B 13 (!), but they also investigated the model extensively with molecular dynamics. In addition they considered the scaled Uhlenbeck-Ford model which has the potential multiplied by an integer. That makes the core more repulsive while keeping the graph contributions still integrable. They also gave a survey of the properties and applications of the UF model.
In this paper we make an extension of the model by adding an attractive part to the potential. One could do this by taking f (r) as a sum of two gaussians Since f (0) = −1 the potential has again a repulsive core. By taking the range a 2 > a 1 the potential gets an attractive tail (for A > 0). By playing with the parameters one can influence the range and the depth of the attractive well in the potential. The price to be paid is considerably more computational effort and each choice of the parameters requires the full evaluation of the virial coefficients. A more severe limitation comes from the fact that we want to have the range a 2 of the attraction not too different from the range a 1 of the repulsion. As the two terms in Eq. (3) have opposite signs, the contributions tend to cancel, making the sum much smaller than the individual terms, which easily leads to numerical errors. This is a delicate problem which generally plagues the evaluation of the virial coefficients, since they are the result of many contributions with uncorrelated signs. A remedy for this danger is to see the difference as a derivative, leading to the Mayer function f (r) = (−1 + A(r/a) 2 ) exp(−(r/a) 2 /2).
This paper is devoted to the evaluation of the virial coefficients resulting from this Mayer function. They become functions B l (A) in the form of a finite polynomial The k-th power of A comes from graphs with at least k edges. The maximum power of A comes from the graph with all m edges occupied with So the expansion in powers of A terminates at the m-th power. For A = 0 only the k = 0 terms contribute and the series Eq. (7) reduces to that for the soft repulsive potential. The model, described by Eq. (6), contains only one free parameter A apart from the scale parameter α. In that sense it is less versatile than the model given by Eq. (5), which contains three extra parameters. But the result as presented in Eq. (7) is more useful. It gives explicitly the A dependence of the virial coefficients B l (A), through the B l,k as a finite series in powers of A.
There is another advantage of the extension Eq. (6). One of the drawbacks of the Uhlenbeck-Ford model is that there is no temperature dependence, like in the hard sphere model. Whereas in the hard sphere model this is a consequence of the fact that the potential is either zero or infinite, the potential has in the Uhlenbeck-Ford model has a distance dependent structure, which cannot be varied by a continuous amplitude. With the choice Eq. (6) the thermodynamic properties become functions of the density n and the parameter A, which mimics the variable 1/T . A large A means a deep well in the reduced potential V (r)/(k B T ) or a low T at fixed V (r).
It is customary in this field to present the results in a dimensionless way, for which we use the second virial coefficient B 2,0 of the repulsive part with d the dimension of the system. With this molecular volume we construct the dimensionless density n * as n * = nv 0 (10) and turn the virial series into with The coefficients B * l,k are dimensionless numbers. They constitute the objects to be calculated in this paper.

Properties of the potential and the parameters
The potential V (r) corresponding to the Mayer function f (r), given by Eq. (6), reads In Fig. 1 we show its behavior for the value A = 1. The potential is zero at the For distances r smaller r 0 the particles feel the repulsive part of the potential and beyond r 0 they are in the attractive well of the potential. The potential V (r) has a minimum where f (r) has a maximum, which occurs at the position r = r m The value of f (r m ) and the depth of the potential V (r m ) are related as In order to illustrate the influence of the parameters n * and A, we consider the configuration where the particles occupy the points of an fcc-lattice. This a closepacked arrangement where the distance of any particle to its nearest neighbors equals the same value b, given by the density n as We compare this distance with the distance r 0 , given by Eq. (14), using the reduced density n * as defined in Eq.
Small A and/or large n * force the fcc distance to be smaller than r 0 and the particles are in the repulsive part of the potential. For larger A or smaller n * they are in the attractive well of the potential. For potentials with a finite hard core radius σ, the system forms a solid before the distance b equals σ. Larger densities n are not possible as they lead to an infinite pressure. Likely this geometric ordering is absent in the Uhlenbeck-Ford model. Also the behavior of the virial coefficients supports this idea as there is no indication of a diverging virial series for a finite radius of convergence.

The Gaussian Integrals
In this Section we illustrate the evaluation of the virial coefficients. In Fig. 2 we show the three graphs contributing to the fourth virial coefficient. The graphs have to be doubly-connected, i.e. they remain connected if one of the vertices is removed [2]. We rewrite the Mayer function Eq. (6) as and set later α = 1. So we have to evaluate the graph contribution as a function of the α ij on the edge (ij). As we shall see this yields a relatively simple polynomial in the α ij . Then perform the differentiations with respect to the α ij and set them equal to α ij = 1. We illustrate the handling of the Gaussian integrals for the fourth virial coefficient. The integrals are of the form where the sum runs over the occupied edges (ij).
In the integral over the vertices we may carry out a permutation of the vertices without changing the contribution. In general such a permutation leads to a different connectivity, i.e. to a different graph. These graphs can be taken together. It amounts to a cancelation of the l! in the denominator of Eq. (2), unless the graph has symmetries. The first graph in Fig The symmetries of a graph are important, not only because of the symmetry number, but also in reducing the computational evaluation. Therefore it is useful to generate the symmetry operations together with the generation of the graph.
The integrand of Eq. (20) is translational invariant, which is exploited by making the shift for i = 0. As result r 0 does not appear in the integrand anymore and the integral over r 0 yields a volume factor in the integration in Eq. (20). By the shift Eq. (21), the exponent is changed into The matrix w i,j is given by (only non-zero i and j).
After the elimination of the integration over r 0 the contribution C 4 becomes The connectivity matrix w i,j is characteristic for the graph. The gaussian integral is evaluated by diagonalization of the matrix w i,j , leading to the eigenvalues λ i on the diagonal. Each eigenmode gives a factor ([2πa 2 ] 3 /λ i ) d/2 with d the dimension of space. Since the product of the eigenvalues equals the determinant we find for C 4 Written out the determinant reads for the graph with k = 4 det(w) = α 01 +α 12 -α 12 0 -α 12 α 12 +α 23 -α 23 0 -α 23 α 03 + α 23 Working out the terms of the determinant yields the function P (α ij ) = α 01 α 12 α 23 + α 12 α 23 α 03 + α 23 α 03 α 01 + α 03 α 01 α 12 .
This function, which is characteristic for the graph, is called the Kirchhoff polynomial. Note that all terms have l − 1 (= 3) factors α ij and no higher powers of the α ij occur. If we draw an edge for every α appearing in a term we get the four spanning trees of the graph. In fact there is a general rule [14] of Kirchhoff linking the graph function to the set of spanning trees of the graph the Kirchhoff polynomial is given by the set of numbered spanning trees of graph.
A spanning tree is a connected graph with the minimal number of edges. So for a graph of l vertices, a spanning tree has l − 1 edges. The full set of spanning trees of 4 vertices are shown in Fig. 3. The 4 terms in Eq (26) correspond to the 4 possible numberings of the spanning tree. The value of the graph function for all α ij = 1 is the number of spanning trees.
Instead of the determinant we use the Kirchhoff polynomial as it contains all the information needed for the differentiations. Spanning trees are easily generated and therefore a convenient way to calculate the Krichhoff polynomial. For the fully occupied graph of l vertices the number of spanning trees is given by the Caley theorem [15] The fully occupied graph has the maximum number of spanning trees. The Kirchhoff polynomial of the fully occupied graph with l vertices, contains also all the information needed for calculation of the other graphs. The Kirchhoff polynomial for graphs with lesser occupied edges follows from the fully occupied graph by setting the α ij = 0 for the empty edges.   Table 1. Virial coefficients B l,0 for pure repulsion.
These numbers agree to all digits with those given in [12]. As they were obtained prior to noticing this paper, it can be seen as a mutual confirmation.

The Lower Derivatives
In order to find the algorithm for the derivatives we first concentrate on the low derivatives as occurring in the first few virial coefficients. The second virial coefficient has a single graph of two vertices connected by an edge. This is also the spanning tree of the graph and the associated polynomial reads For the virial coefficient we find Thus we have the two coefficients, B 2,0 given by Eq. (9) and The third virial coefficient also involves a single graph, the fully occupied graph with three vertices connected by three edges. Numbering the edges 1,2 and 3 we obtain the graph polynomial P (α 1 , α 2 , α 3 ) = α 1 α 2 + α 2 α 3 + α 3 α 1 .
The third virial coefficient is then found as Performing the differentiations and setting the α's equal to 1 yield for B * from which the components B * 3,k follow as the coefficients of the powers in A. Generally, let κ, λ, µ be a choice of three edges out of the k edges of the graph. Then we have the relations This sequences of relations is general and holds for all graphs. Each higher derivative follows from the previous one. The number of terms grows since we have to calculate the derivatives of a (negative) power of P . Each term of the lower derivative yields a number of terms in precisely the same way as the partitions of a set with q elements follow from those with q − 1 elements. So the systematics from the sequence is clear: consider for the last Eq. (33) all the partitions of a set κ, λ, µ into groups and take the products of these groups of derivatives. The dimensional weights of these terms depend only on the number of groups (or the number of factors). The sequence also shows that the number of terms grows very rapidly with the number of derivatives involved, as rapid as the number of partitions of a set grows with the number of elements in the set.
The partitions are independent of the polynomials P and they can be generated once for all graphs and virial coefficients. The problem is that the number increases so fast that memory problems arise sooner or later for the higher number of edges. Therefore we have to design strategies for handling this problem.

Algorithms for the higher derivatives
We outline here two complementary methods of evaluating the coefficients B l,k . The first uses the partitions and the second computes the derivatives recursively. All numerical calculations were carried out for dimension d = 3.

The method using partitions
Here the basic ingredients are the derivatives of the Kirchhoff polynomial P . Their values are integers: the number of terms that survive after the differentiation. Since all q-th order derivatives vanish for q >= l, there are not so many non-vanishing and to begin with, they can be listed for each graph.
Then we have to carry out three summations: 1. The basic summation over the graphs qualifying for the virial coefficient.
2. For the coefficient B l,k all the choices of k edges out of the set of edges of the graph.
3. For each choice of k edges the summation over all partitions of the k edges. Each partition contributes a product of the listed derivatives of P . This is straigthforward and fast, but gets time comsuming and inaccurate for the highest derivatives. The errors occur because the partitions get a sign depending on the number of factors: even numbers of factors have a positive weight and odd numbers a negative weight, while their contributions are often comparable. These errors are avoided by integer arithmetic. Time and accuracy is gained by using the symmetries of the graph. Several choices of k edges can be equivalent due to the symmetry of the graph. E.g. any choice of two edges in the case of the third virial coefficient gives the same result. Calculating one of them and multiplying with the number of equivalent choices is sufficient.
Without integer arithmetic the method works satisfactorily up to 11th derivative, with integer arithmetic we could extend that by a few more.
In some cases the sum over the various derivatives simplifies due to the following rule: Here σ is a set of indices (κ, λ, · · ·) with q terms. The proof of this rule is based on the fact that one can select in l − 1 over q ways an edge in each term of the spanning trees. All these selections contribute 1 (after setting the α i = 1) and the total involves the number of terms P in the spanning tree of the graph. Doing the summation over the derivatives in Eqns. (33) simplifies the outcome for to This rule can be used for the first term on the right hand sides of the Eqns. (33). For the first power of A one finds The proof of this relation follows from the property Eq. (38).

The recursive method
The second method computes the derivatives of P −d/2 recursively. Let again σ be a choice κ, λ, · · · , µ, ν of k edges of the graph. As the order of the differentiations does not matter we order them in increasing such that ν is larger than the preceding ones. We can compute the derivative with respect to the last chosen edge ν from the derivative of the predecessing choice κ, λ, · · · , µ. These derivatives are ratios of two polynomials and can be written as The denominator is explicit and the numerator obeys the recursion To make the recursion complete we set R 0 = 1.
To illustrate the recursion we take the case l = 3 as example with P given by Eq. (31). The first level then gives The next level leads to Here we encounter a new element. After a differentiation is carried out, we may set the corresponding α = 1. So the factor P in R 2 reduces to P = 1 + 2α 3 and R 2 becomes a function of the last remaining variable α 3 .
The final R 3 is a number Note that this result agrees with the values given in Eq. (33). Characteristic for this method is that the polynomials increase in length with the number of differentiation k due to the multiplications, but they shrink because the increasing number of α's which can be put equal to 1. The latter property is strengthend due to our choice to carry out the differentiations in the order of increasing edge index. That means that all the α's with index lower than the differentiated α may be put equal to 1, since they will not occur in a further choice. When all the k differentiations are carried out the resulting polynomial is a number.
While the case, where three differentiations have to be carried out, can be done by hand, higher orders have to be programmed. The problem is to add and multiply polynomials. Therefore the terms in the polynomial have to be coded. They have a coefficient in front and a row of powers of the variables α. The variables α which are set equal to 1 get a power 0. Multiplying two terms implies multiplication of the coefficients and addition of of the exponents. Since the factors involving P have at most a power 1, in each multiplication the maximum power is raised by 1.
It means that powers never exceeds k. (E.g. R 2 in Eq. (44) has as highest power 2.) So powers form a row of k digits. The digits are limited by the order k of the differentiation. The factors P and ∂ ν P contain the α i to atmost order 1. Therefore the maximum power in R increases in each recursion by 1. After setting α ν = 1 the number of remaining α i decreases by 1. So the expression for R k first increases in complexity by increase of k and finally reduces to a single number when k reaches the maximum value m of edges. In symmary the highest power in R k is therefore k and the number of remaining variables is m − k. This makes the method efficient for the larger k. A graph with q edges contributes to the B l,k with k ≤ q.
The two methods are complementary: the first spends the most time in calculating the highest derivatives, the second is slowest in the medium derivatives and speads up towards the largest derivatives. The two have to give the same result and this is a superb check on the calculation. For l = 6 the polynomial method spends too much time on the last few graphs with 14 and 15 edges to complete the calculation, while the method using the partitions just barely makes it in an acceptable running time.
Below we have listed the coefficients for l up to 6 and k up to 10 for dimensions d = 3.   Table 2. Polynomial coefficients B l,k for the extended UF-model. The 5 further coefficients B * 6,k for k = 11, 12, 13, 14 and k = 15 are respectively: -16.65300, -2.69820, -0.59634, 0.14819 and -0.04533. The last collumn refers to the approximation described in the Appendix.
In order to see what these polynomial coefficients mean for the virial coefficient B * l (A) we plot these coefficiens as function of A in Fig. (4). Apart from B * 2 (A), which depends linearly on A, they drop off rapidly to small values in the range A = 0.2 to A = 0.6. Beyond A = 0.6 they evolve to larger negative values. The range A = 0.2 to A = 0.6 of small virial coefficients point to a delicate interplay of the terms is the polynomial expression for B * l (A). So the B * l,k have to be computed with high precision. In order to see how small the values become near A=0. 5 Table 2, gives the values of B * l,k following from the approximation described in the Appendix.

Phase Diagram
Speculating on the phase diagram of the extended UF model, we note that the potential does not have a hard core with a finite range or a sharp deep well, which may induce the geometric ordering in a crystal. So we do not expect that the model has a transition to an ordered solid phase at high densities.
On the other hand the attraction may lead to a gas-fluid phase transition at intermediate densities. As a signal of this transition one would see the appearance of a van der Waals loop in the pressure. Such a loop is the result of an interplay of a negative second virial coefficient bending the pressure down and positive higher virial coefficients which turn the pressure upwards again at higher densities.
We note that in Table 2 the B l,m , with m = l(l − 1)/2 the highest power of A, all are negative. That means that for large A the virial coefficients (beyond the ideal gas term) become negative and that there is no stability in the pressure at high pressures. At intermediate values of A, 0.4 < A < 0.6, the highest B l,m is positive for l = 3 and l = 6. For these cases a van der Waals loop may occur.
The onset of the loop, the critical point, is found from the conditions Using the virial series for the pressure, the equations get the form Using Eq. (29) for B * 2 (A) and Eq. (33) for B * 3 (A) gives a cubic equation for A, with the solution A c = 0.43647. For the virial series terminating at l = 6 one has to solve the equation numerically, with result A c = 0.42591 In Fig. (5) we show a few "isothermes", the critical and two subcritical exhibiting a vander Waals loop, with the equal area construction for the coexisting phases. This gives an indication of the gas-fluid transition and its location in the phase space. Of course one cannot derive the critical singularities from a finite virial series. It is likely that the critical point of the model with attraction is in the class of the 3d-Ising model.

A Connection to standard potentials
In order to get an impression of the length scale a and the attraction amplitude A, we compare the a few virial coefficients of the model, with those of 2n − n Lennard-Jones potentials. For n = 6 this potential is often used for simple molecules like the noble gases. They are of the form [16] is the depth of the potential well and σ the range of the interaction. can be used for scaling the temperature T to the dimensionless T * as Let f (r) be the associated Mayer function according to Eq. (3), then the second virial coefficient B 2 can be written as B(T * ) is a dimensionless function of T * . Comparing this with the expression Eq. (37) for the second virial coefficient we see that one can tune the parameter a/σ such that the two values are equal with one proviso: the parameter A has to be chosen such that both virial coefficients have the same sign. That means that for higher T * , where the Lennard-Jones second virial coefficient is positive, A < 1/3 and that for the lower temperatures A > 1/3. With this restriction A is a still a free parameter which can be tuned such that also the third coefficients of the representations coincide.
The third virial coefficient is given by Using the fourier transform the third virial coefficient reads The ratio is independent of the scale σ and may be used for matching the third virial coefficient. The corresponding ratio for the soft potential is contained in Eq. (31) and reads explicitly C(T * ) This gives a relation between T * and A. Note that this relation does not involve the ranges of the two models which are compared. In Fig. (6) we show for a few choices of n the value of A that corresponds to T * . This relation implies an estimate for the critical temperature T * c (n). We find T * c (8) = 0.95, T * c (7) = 1.16, T * c (6) = 1.54, T * c (5) = 2.36, T * c (4) = 4.56, which are reasonable estimates, given the crudeness of the match and which show the proper trend. The vertical line corresponds to the critical value A c = 0.4259

Conclusion
We have extended the Uhlenbeck-Ford model with an attractive part, having an amplitude A. The corresponding potential is drawn for A = 1 in Fig. 1 m = l(l − 1)/2. The calculation of the coefficients of these power series is considerably more involved than that of B l for A = 0 (pure repulsion). We have determined B l (A) for l up to 6.
The key quantity in the calculation is the Kirchhoff polynomial P (α i ), defined in Section 3. It is a polynomial in the variables α i , based on the spanning trees of the graph. Each of the l − 1 edges of a spanning tree contributes a factor α i if the i-th edge is present in the graph. For B l (A) one needs the derivatives of P −d/2 . We have employed two methods: one using the partitions of the α i and one using recursively the derivates of the function P (α i ). The mutual agreement of these independent methods is a guarantee of the correctness of the coefficients.
One may expect that with stronger attraction the model forms a fluid phase. Indeed for intermediate A we observe the formation of a van der Waals loop with a critical value A c = 0.246. An intruiging question is the behaviour for large values of A of the functions B l (A). The sofar calculated highest coefficients B l,m become negative, which implies that the functions B l (A) becomes negative for large A. That leads to an instability in the virial series through negative values of the pressure.
We have made a link between the amplitude A of the attraction and 1/T * in Lennard-Jones type interaction.
An approximation scheme has been developed (see the Appendix) which works well for the lower derivatives.
We close with the remark that the graph contributions to the pair correlation function are also exactly calculable as a power series in the density, but this is another project, more complicated than the present one.

A Approximation
As mentioned the number of partitions of a set increases very rapidly with the number of elements (edges of the graph). Therefore the computation of the virial coefficients spends most of the time calculating the contribution of the (almost) fully occupied graphs. We also observed that the various derivatives of a set of edges vary little with the specific edges in the set. Only the size q of the set is the main ingredient. Note that for the average value one can use the rule Eq. (38).
Based on this idea we have designed an fast scheme, which give an indication of the magnitude of the contributing terms. Replacing the individual derivatives by the average over a set of q edges, many partitions give the same value for the derivative. Thus we can lump these partitions together. In fact all the permutations of the edges which keep the sizes of the bins equal, yield the same contribution. Their number is the multinomial N ({σ}) = σ! (σ 1 !σ 2 ! · · ·)(p 1 !p 2 ! · · ·) .
Here {σ} = σ 1 , σ 2 · · · is a set of bin sizes with σ 1 + σ 2 · · · = σ and p 1 the number of bins of size 1, p 2 of size 2 · · ·, with p 1 + 2p 2 + · · · = σ. The contribution of the selection of partitions then equals This is an enormous reduction in the number of partitions. For instance for k = 15 edges to be distributed over 4 bins, one has 42355950 detailed partitions which all lead to the same average value. In fact one needs only the partitions of the edges over bins, no matter which edge is in which bin (the partitions of a number). If one wants to do better than this lowest approximation, one must take the fluctuations into account. This can be done by giving each bin of edges the sum of the average value of the derivative and the deviation thereof (the fluctuation). The lowest approximation is taking only average values. The first approximation is taking one fluctuating bin combined with averages of the others (yields a vanishing contribution), the second by taking two fluctuating bins combined with averages etc. So one still needs the full partitions of the fluctuating bins, but that is limited as long as their number is small. This gives a quick estimate of the value of B * l,k provided that k is not too large. In fact if one takes into account the products of fluctuation to order n the first cofficients B * l,k with k ≤ n becomes exact. In Table 2, last column, we have given the approximate values, taking the fluctuations into account to 4th order. For all other values in the Table the approximation gives virtually the exact value. Beyond k = 10 and l = 6 the approximation gets inaccurate, but these coefficients are less important for A < 1.
As Table 2 shows the approximated B * 6,k are quite close to the exact values. But these small differences can build up to important differences in the value of B * l (A), as the individual terms B * 6,k A k are much larger than the sum B * l (A) for A 0.5, which is is an interesting value of A (see next Section (6).