text:
# Dimers on a simple-quartic net with a vacancy

## 1 Introduction

The problem of enumerating close-packed dimers on a finite simple-quartic net was solved by Temperley and Fisher  and by Kasteleyn  in 1961. An outstanding related but yet unsolved problem is the enumeration of dimers on a net with vacant sites . Here, we consider this vacancy problem when a single vacant site occurs on the boundary.

The difficulty associated with the vacancy problem is that, while the determinant whose square root yieids the dimer generating function can be written down using the Kasteleyn formulation , its evaluation is difficult. In 1974 Temperley  reported an intriguing bijection relating close-packed dimer coverings with spanning tree configurations on two related lattices. This offers an alternate approach to the vacancy problem since spanning trees can be enumerated by standard means. The Temperley bijection has been of renewed recent interest and extended to graphs with certain weighted and/or directional edges . Here, we use an extension of the Temperley bijection due to Kenyon, Propp, and Wilson  to study the vacancy dimer problem.

Our first result is that, for a simple-quartic net with free boundaries and one fixed vacant site located at certain specific sites on the boundary, the dimer generating function is independent of the position of the vacancy. The exact generating function for close-packed dimers on this net is then deduced from that of the spanning trees. In view of the connection with the conformal field theory  and current interests in finite-size analyses for two-dimensional lattice models , we next carry out finite-size analyses for the weighted spanning tree solution as well as that of the FKT solution. It is found that a logarithmic correction term arises in the large-size expansion in the case of the vacancy problem with free boundaries, a term which is absent in the expansion of the FKT solution. A concrete example demonstrating this difference is given. We also find that the occurrence of a vacancy yields a new central charge $`c=2`$ in the language of the conformal field theory.

The organization of this paper is as follows: To make the paper self-contained, we restate and establish in section 2 the bijection due to Temperley  between spanning tree and dimer configurations, as well as an extended version of the bijection due to Kenyon, Propp, and Wilson . This extended Temperley bijection permits us to establish in section 3 the independence of the dimer generating function on the location of the vacancy when it occurs at specific boundary sites. The explicit expression of the generating function for a simple-quartic net is then obtained. In section 4 we carry out finite-size analyses for weighted spanning trees and the FKT dimer solution. Specializing the results to dimer enumerations, we find a logarithmic correction term which is unique to the vacancy problem with free boundaries. Discussions and a summary of our findings are given in section 5.

## 2 The extended Temperley bijection

For definiteness we consider a simple-quartic net (lattice) with free boundaries, although much of the results of this section also hold for more general planar graphs .

First we restate the Temperley bijection. Starting from a $`L_1\times L_2`$ simple-quartic net $`G`$ with free boundaries, one constructs a dimer lattice $`G_D`$ by i) adding a new site at the midpoint of each edge of $`G`$, ii) inserting in each internal face of $`G`$ a new site connected to the midpoints of the 4 edges of $`G`$ surrounding it, and iii) removing one corner site of the resulting lattice and its incident edges on $`G_D`$. Thus, $`G_D`$ has a total of $`(2L_11)(2L_21)1`$ sites consisting of the original $`L_1L_21`$ sites of $`G`$, which we call the odd sites, and the remaining $`(2L_11)(2L_21)L_1L_2`$ new sites, which we call the even sites. An example of this construction for $`L_1=L_2=3`$ is shown in Figs. 1(a) and 1(b).

A spanning tree is a collection of connected edges of $`G`$ which does not form closed circuits and covers all sites. Then we have the

Temperley bijection: There exists a one-one correspondence between spanning tree configurations on $`G`$ and dimer configurations on $`G_D`$.

To see that the bijection holds, one observes that to each spanning tree configuration on $`G`$, one can construct a unique dimer configuration on $`G_D`$ by first laying a dimer along each tree edge, starting from the edge(s) covering the corner site of $`G_D`$ which has (have) been removed, and proceed along the spanning tree edges in an obvious fashion. After laying dimers along all tree edges, the remaining sites of $`G_D`$ can then be covered by dimers in a unique way . Conversely, starting from each dimer configuration on $`G_D`$, one constructs a unique tree configuration on $`G`$ by drawing bonds (tree edges) along dimers originating from all odd sites. These bonds cannot form close circuits, since otherwise they would have enclosed an odd number of sites of $`G_D`$ which is not permitted in close-packed dimer configurations. This process leads to a unique tree configuration on $`G`$. This completes the proof. An example of the Temperley bijection is shown in Figs.2(a) and 2(b).

Kenyon, Propp, and Wilson  have shown that the Temperley bijection holds more generally for graphs with certain weighted and/or directed edges. For our purposes, however, we shall confine ourselves to the original Temperley bijection as stated in the above, with the step iii)replaced by one of removing an odd site on the boundary together with its incident edges. An example of constructing such a $`G_D`$ is shown in Fig.1(c). The proof of the bijection between tree configurations on $`G`$ and dimer coverings on $`G_D`$goes through as before, and we are led to the

Extended Temperley bijection (Temperley-Kenyon-Propp-Wilson): There existsa one-one correspondence between spanning trees on $`G`$ and dimer coverings on any $`G_D`$ constructed from $`G`$ by removing any boundary odd site and its incident edges in step iii) of the construction described in the above.

An example of such a bijection is shown in Figs.2(a) and 2(c).

Remark: The extended Temperley bijection does not hold for dimer lattices $`G_D`$ containing an interior vacancy. In that case while each spanning tree can still be mapped into a unique dimer configuration as before, there exist dimer configurations which cannot be mapped into spanning trees. These are dimer coverings with no dimers laying on any of the 4 edges of $`G`$ incident to the defect site.

## 3 Dimer lattice with a vacant boundary site

### 3.1 Dimer generating function

The dimer generating function for $`G_D`$ is

$$Z(G_D;x_1,x_2)=\underset{\mathrm{dimer}\mathrm{config}.}{}x_1^{n_1}x_2^{n_2}$$
(1)
where the summation is taken over all dimer covering configurations, $`x_1`$ and $`x_2`$ are, respectively, the weights of horizontal and vertical dimers., and $`n_1`$ and $`n_2`$ are, respectively, the number of horizontal and vertical dimers. Clearly, we have

$$Z(G_D;1,1)=\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{dimer}\mathrm{configurations}\mathrm{on}G_D.$$
(2)
Consider two different dimer lattices $`G_D`$ and $`G_D^{}`$ obtained from $`G`$ as described in the above, namely, by removing different boundary odd sites. Then we have the following equivalence:

Proposition 3.1.1:

$$Z(G_D;x_1,x_2)=Z(G_D^{};x_1,x_2)$$
(3)
for any $`G_D`$ and $`G_D^{}`$.

Proof: The extended Temperley bijection dictates that there is a one-one correspondence between spanning tree configurations on $`G`$ and dimer configurations on any $`G_D`$. It follows that there is a bijection between dimer coverings on $`G_D`$ and $`G_D^{}`$, and that the summation in (1) on $`G_D`$ can be considered as taken over all spanning tree configurations on $`G`$.

For each spanning tree configuration $`T`$ of $`G`$, the dimer weight in the summand in (1) consists of two factors,

$$x_1^{n_1}x_2^{n_2}=W_o(T;x_1,x_2)W_e(T;x_1,x_2),$$
(4)
where $`W_o`$ is the product of the weights of those dimers originating from odd sites, and $`W_e`$ is the product of the weights of those dimers covering two even sites. For the two dimer coverings of $`G_D`$ and $`G_D^{}`$ corresponding to the same$`T`$, their $`W_e`$ factors are the same by definition. Their $`W_o`$ factors are also the same since, even though the respective dimer positions may be shifted, they lay along the same spanning tree edges hence carry the same weights. It follows that summations on the lhs and rhs of (3) are identical term by term, and the proposition is proved. $``$

Remark: Proposition 3.1.1 holds more generally for arbitrary planar $`G`$ and its related $`G_D`$. Since the overall dimer weight can always be factorized into the product $`W_oW_e`$ as in (4), the proof of the proposition goes through as presented.

We next consider the generating function of weighted spanning trees for the $`L_1\times L_2`$ simple-quartic $`G`$. Assign weights $`x_1`$ and $`x_2`$, respectively, to edges in the horizontal and vertical direction. Then, the weighted spanning tree generating function is

$$T(G;x_1,x_2)=\underset{T}{}x_1^{n_1}x_2^{n_2}$$
(5)
where the summation is taken over all spanning tree configurations $`T`$ on $`G`$ and, as in (1), $`n_1`$ and $`n_2`$ are the numbers of edges in the spanning tree in the respective directions. Particularly, we have

$$T(G;1,1)=\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{spanning}\mathrm{tree}\mathrm{configurations}\mathrm{on}G.$$
(6)
From the extended Temperley bijection, it is clear that we can also write (5) as

$$T(G;x_1,x_2)=\underset{T}{}W_o(T;x_1,x_2)$$
(7)
where $`W_o`$ is the factor in (4) for the dimer covering on any $`G_D`$. It is seen from (7) that if all dimers covering even sites of $`G_D`$ have the weight $`1`$, namely $`W_e=1`$, then the dimer generating function is simply $`T(G;x_1,x_2)`$.

More generally for a simple-quartic $`G`$ of size $`L_1\times L_2`$ and the related $`G_D`$, we have the equivalence:

Proposition 3.1.2:

$$Z(G_D;x_1,x_2)=x_1^{L_1(L_21)}x_2^{L_2(L_11)}T(G;\frac{x_1}{x_2},\frac{x_2}{x_1}).$$
(8)
Proof: From the construction of $`G_D`$ we note that there are a total of $`L_2(L_11)`$ even sites located at midpoints of horizontal edges of $`G`$. We call these the $`H`$ sites. Similarly, there are $`L_1(L_21)`$ even sites of $`G_D`$ located at midpoints of vertical edges of $`G`$, which we call the $`V`$ sites.

The $`H`$ and $`V`$ sites can be covered by either horizontal or vertical dimers. Let $`N_{Hh}`$ ($`N_{Hv}`$) be the number of horizontal (vertical) dimers covering the $`H`$ sites. Then we have

$$N_{Hh}+N_{Hv}=L_2(L_11).$$
(9)
Likewise, we have

$$N_{Vh}+N_{Vv}=L_1(L_21),$$
(10)
where $`N_{Vh}`$ ($`N_{Vv}`$) is the number of horizontal (vertical)dimers covering the $`V`$ sites. In these notation, we can rewrite the spanning tree and dimer generating functions as

$`T(G;x,y)`$ $`=`$ $`{\displaystyle \underset{T}{}}x^{N_{Hh}}y^{N_{Vv}}`$
$`Z(G_D;x_1,x_2)`$ $`=`$ $`{\displaystyle \underset{T}{}}x_{1}^{}{}_{}{}^{N_{Hh}+N_{Vh}}x_{2}^{}{}_{}{}^{N_{Hv}+N_{Vv}},`$ (11)
where we have used the one-one correspondence between spanning tree configurations $`T`$ and dimer configurations. The Proposition 3.1.2 now follows after eliminating $`N_{Vh}`$ and $`N_{Hv}`$ in (11) using (9) and (10). $``$

### 3.2 Dimer enumerations

For a simple-quartic $`G`$ of size $`L_1\times L_2`$ with free boundaries, the generating function (7) for weighted spanning trees has been evaluated  and is given by

$`T(G;x_1,x_2)={\displaystyle \frac{1}{L_1L_2}}{\displaystyle \underset{m=0}{\overset{L_11}{}}}{\displaystyle \underset{n=0}{\overset{L_21}{}}}\left[2x_1\left(1\mathrm{cos}{\displaystyle \frac{m\pi }{L_1}}\right)+2x_2\left(1\mathrm{cos}{\displaystyle \frac{n\pi }{L_2}}\right)\right],`$
$`(m,n)(0,0).`$ (12)
Now the dimer lattice $`G_D`$ is of size $`M\times N`$ with a boundary vacancy, where

$$M=2L_11,N=2L_21.$$
(13)
In terms of the dimer lattice sizes $`M,N`$, we thus have, after using Proposition 3.1.2 and (12),

$`Z_{\{M\times N1\}}(x_1,x_2)=x_1^{(M1)/2}x_2^{(N1)/2}`$
$`\times {\displaystyle \underset{m=1}{\overset{\frac{M1}{2}}{}}}{\displaystyle \underset{n=1}{\overset{\frac{N1}{2}}{}}}[4x_{1}^{}{}_{}{}^{2}\mathrm{cos}^2{\displaystyle \frac{m\pi }{M+1}}+4x_{2}^{}{}_{}{}^{2}\mathrm{cos}^2{\displaystyle \frac{n\pi }{N+1}}].`$ (14)
Here, we must have $`MN=`$ odd to admit dimer coverings. The subscript $`\{M\times N1\}`$ in (14) reminds us that the enumeration is for an $`M\times N`$ net with one boundary odd site removed. Note that the factor $`m=n=0`$ is excluded in the product in (14). This expression is to be compared with the enumeration of dimers on an $`M\times N`$ simple-quartic net without vacancies. For $`M`$ and $`N`$ both even, for example, the expression is

$$Z_{\{M,N\}}(x_1,x_2)=\underset{m=1}{\overset{M/2}{}}\underset{n=1}{\overset{N/2}{}}\left[4x_1^2\mathrm{cos}^2\left(\frac{m\pi }{M+1}\right)+4x_2^2\mathrm{cos}^2\left(\frac{n\pi }{N+1}\right)\right].$$
(15)
## 4 Finite-size analyses

Finite-size expansions of physical quantities associated with two-dimensional lattice models have been of current interest both in physics  and in mathematics . For the dimer problem Kenyon  has recently deduced very general results on the leading terms of the asymptotic expansion of the dimer enumeration for rectilinear lattices with free boundaries of any shape. Alternately, one can obtain expansions for regular lattices, in principle to all orders, by analyzing known exact expressions. This is the approach we now use.

### 4.1 Spanning tree generating function

Consider first the generating function (12) for the fully weighted spanning trees. For large $`L_1`$ and $`L_2`$, we expect to have

$$\frac{1}{L_1L_2}\mathrm{ln}T(G;x_1,x_2)=f_{\mathrm{bulk}}(x_1,x_2)+f_c(x_1,x_2)$$
(16)
where $`f_{\mathrm{bulk}}`$ is the per-site bulk free energy and $`f_c`$ is the correction containing terms of the order of $`L_1^1`$, $`L_2^1`$ and higher. Using (12), we find the bulk free energy

$`f_{\mathrm{bulk}}(x_1,x_2)`$ $``$ $`\underset{L_1,L_2\mathrm{}}{lim}{\displaystyle \frac{1}{L_1L_2}}\mathrm{ln}T(G;x_1,x_2)`$ (17)
$`=`$ $`{\displaystyle \frac{1}{\pi ^2}}{\displaystyle _0^\pi }𝑑\theta {\displaystyle _0^\pi }𝑑\varphi \mathrm{ln}\left[2x_1(1\mathrm{cos}\theta )+2x_2(1\mathrm{cos}\varphi )\right]`$
$`=`$ $`{\displaystyle \frac{4}{\pi }}{\displaystyle _0^{\pi /2}}𝑑\varphi \mathrm{ln}\left(\sqrt{x_1+x_2\mathrm{sin}^2\varphi }+\sqrt{x_2}\mathrm{sin}\varphi \right),`$
where the last line is obtained by carrying out the $`\theta `$ integration. The computation of correction terms $`f_c`$ for products of the form of (12) is standard . Particularly, one has

$$f_{\mathrm{bulk}}(1,1)=\frac{4}{\pi }G$$
(18)
where $`G`$ is the Catalan constant given by

$$G=13^2+5^27^2+\mathrm{}=\mathrm{0.915\; 965\; 594}\mathrm{}.$$
(19)
To compute (16) we proceed as follows. Take out a factor $`x_1`$ from each of the $`L_1L_21`$ factors in(12) and split the product into 3 parts to take care of the exclusion of the $`m=n=0`$ factor.

We have

$$T(G;x_1,x_2)=\left(L_1L_2\right)^1x_1^{L_1L_21}\left(T_0T_1T_2\right)$$
(20)
where

$$T_0=\underset{m=1}{\overset{L_11}{}}\underset{n=1}{\overset{L_21}{}}F(m,n),T_1=\underset{m=1}{\overset{L_11}{}}F(m,0),T_2=\underset{n=1}{\overset{L_21}{}}F(0,n)$$
(21)
with

$$F(m,n)=2\left(1\mathrm{cos}\frac{m\pi }{L_1}\right)+2\tau \left(1\mathrm{cos}\frac{n\pi }{L_2}\right)$$
(22)
and

$$\tau =x_2/x_1.$$
(23)
Using the identity

$$\underset{n=1}{\overset{N1}{}}\left[2\mathrm{cosh}2\theta 2\mathrm{cos}\frac{n\pi }{N}\right]=\frac{\mathrm{sinh}2N\theta }{\mathrm{sinh}2\theta },$$
(24)
and its $`\theta 0`$ limit,

$$\underset{n=1}{\overset{N1}{}}\left[22\mathrm{cos}\frac{n\pi }{N}\right]=N,$$
(25)
we find

$$T_0=\underset{n=1}{\overset{L_21}{}}\frac{\mathrm{sinh}(2L_1\theta _n)}{\mathrm{sinh}2\theta _n},T_1=L_1,T_2=\tau ^{L_21}L_2,$$
(26)
with $`\theta _n`$ given by

$$\mathrm{cosh}2\theta _n=1+\tau \left(1\mathrm{cos}\frac{n\pi }{L_2}\right),$$
(27)
or, explicitly,

$`\theta _n=F\left({\displaystyle \frac{n\pi }{2L_2}}\right)`$ $``$ $`\left[\mathrm{cosh}^1\left(1+2a_n^2\right)\right]/2`$ (28)
$`=`$ $`\mathrm{sinh}^1a_n`$
$`=`$ $`\mathrm{ln}\left(a_n+\sqrt{1+a_n^2}\right),`$
where

$$a_n=\sqrt{\tau }\mathrm{sin}\frac{n\pi }{2L_2}.$$
(29)
Substituting (26) into (12), we thus obtain

$$T(G;x_1,x_2)=x_1^{L_1L_21}\tau ^{L_21}\underset{n=1}{\overset{L_21}{}}\frac{\mathrm{sinh}(2L_1\theta _n)}{\mathrm{sinh}2\theta _n}.$$
(30)
The product in the denominator in (30) can again be evaluated using (24) as

$`{\displaystyle \underset{n=1}{\overset{L_21}{}}}\mathrm{sinh}^22\theta _n`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{L_21}{}}}(\mathrm{cosh}2\theta _n1){\displaystyle \underset{n=1}{\overset{L_21}{}}}(\mathrm{cosh}2\theta _n+1)`$ (31)
$`=`$ $`{\displaystyle \underset{n=1}{\overset{L_21}{}}}\left[\tau \left(1\mathrm{cos}{\displaystyle \frac{n\pi }{L_2}}\right)\right]{\displaystyle \underset{n=1}{\overset{L_21}{}}}\left[2+\tau \tau \mathrm{cos}{\displaystyle \frac{n\pi }{L_2}}\right]`$
$`=`$ $`L_2\left({\displaystyle \frac{\tau }{2}}\right)^{2(L_21)}\left({\displaystyle \frac{\mathrm{sinh}2L_2\alpha }{\mathrm{sinh}2\alpha }}\right),`$
where $`\alpha `$ is given by

$$\mathrm{cosh}2\alpha =1+2\tau ^1$$
(32)
or $`\mathrm{sinh}\alpha =1/\sqrt{\tau }`$, or explicitly,

$$\alpha =\mathrm{ln}\left(\sqrt{\tau ^1}+\sqrt{1+\tau ^1}\right).$$
(33)
Combining these results, we obtain from (30) the expression

$$T(G;x_1,x_2)=x_1^{L_1L_21}\left(\frac{\mathrm{sinh}2\alpha }{L_2\mathrm{sinh}2L_2\alpha }\right)^{1/2}\underset{n=1}{\overset{L_21}{}}\left[2\mathrm{sinh}\left(2L_1\theta _n\right)\right].$$
(34)
Taking the logarithm, we obtain

$`\mathrm{ln}T(G;x_1,x_2)=(L_1L_21)\mathrm{ln}x_1+2L_1{\displaystyle \underset{n=1}{\overset{L_21}{}}}\theta _n+{\displaystyle \underset{n=1}{\overset{L_21}{}}}\mathrm{ln}\left(1e^{4L_1\theta _n}\right)`$
$`L_2\alpha {\displaystyle \frac{1}{2}}\mathrm{ln}\left(1e^{4L_2\alpha }\right){\displaystyle \frac{1}{2}}\mathrm{ln}L_2+{\displaystyle \frac{1}{2}}\mathrm{ln}(\mathrm{sinh}2\alpha ).`$ (35)
For large $`L_1`$ and $`L_2`$ with the ratio $`L_1/L_2`$ finite, the first two terms in (35) contribute to the bulk free energy $`f_{\mathrm{bulk}}(x_1,x_2)`$ given in (17). To carry out the summations in (35), we use the Euler-MacLaurin summation formula given by

$`{\displaystyle \underset{n=1}{\overset{N}{}}}f(a+n\delta )`$ $`=`$ $`{\displaystyle \frac{1}{\delta }}{\displaystyle _a^{a+N\delta }}f(x)𝑑x+{\displaystyle \frac{1}{2}}\left[f(a+N\delta )f(a)\right]`$ (36)
$`+{\displaystyle \frac{\delta }{12}}\left[f^{}(a+N\delta )f^{}(a)\right]+O(\delta ^3).`$
With $`a=0`$, $`N=L_2`$, $`\delta =\pi /2L_2`$, and $`f(x)=F(x)`$ defined in (28), one has

$`{\displaystyle \underset{n=1}{\overset{L_21}{}}}\theta _n`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{L_2}{}}}\theta _n\theta _{L_2}`$ (37)
$`=`$ $`{\displaystyle \frac{L_2}{2}}\left[f_{\mathrm{bulk}}(x_1,x_2)\mathrm{ln}x_1\right]{\displaystyle \frac{1}{2}}\mathrm{ln}(\sqrt{1+\tau }+\sqrt{\tau })`$
$`{\displaystyle \frac{\pi \sqrt{\tau }}{24L_2}}+O(L_2^3).`$
For the second summation in (35), we follow the manipulation in  to write

$$\underset{n=1}{\overset{L_21}{}}\left[1e^{4L_1\theta _n}\right]=\underset{n=1}{\overset{\mathrm{}}{}}\left(1e^{2nL_1\pi \sqrt{\tau }/L_2}\right)+O(L_2^{2+ϵ})$$
(38)
for some $`0<ϵ<2`$, with $`ϵ0`$ when $`L_2\mathrm{}`$. Putting the results together, we find the finite-size correction

$$f_c(x_1,x_2)=\frac{c_1(x_1,x_2)}{L_1}+\frac{c_2(x_1,x_2)}{L_2}+\frac{c_3(x_1,x_2)}{L_1L_2}+o\left(\frac{1}{L_1L_2}\right),$$
(39)
where

$`c_1(x_1,x_2)`$ $`=`$ $`\mathrm{ln}\left(\sqrt{\tau ^1}+\sqrt{1+\tau ^1}\right)`$
$`c_2(x_1,x_2)`$ $`=`$ $`\mathrm{ln}\left(\sqrt{\tau }+\sqrt{1+\tau }\right)`$
$`c_3(x_1,x_2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}L_2+\mathrm{ln}2+{\displaystyle \frac{1}{4}}\mathrm{ln}{\displaystyle \frac{(x_1+x_2)}{x_1^3x_2^2}}{\displaystyle \frac{\pi \sqrt{\tau }L_1}{12L_2}}`$ (40)
$`+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{ln}\left(1e^{2n\pi \sqrt{\tau }L_1/L_2}\right).`$
Particularly, for $`x_1=x_2=1`$, the expression $`f_c(1,1)`$ given by (39) reduces to the one given in  and .

Despite its appearance, the expression for $`c_3(x_1,x_2)`$ is actually symmetric in $`\{x_1,L_1\}\{x_2,L_2\}`$, a fact can be seen from the identity

$`{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\mathrm{ln}\left(1e^{2m\pi L_2/\sqrt{\tau }L_1}\right){\displaystyle \frac{1}{2}}\mathrm{ln}L_1{\displaystyle \frac{1}{4}}\mathrm{ln}x_2{\displaystyle \frac{\pi L_2}{12\sqrt{\tau }L_1}}`$ (41)
$`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{ln}\left(1e^{2n\pi \sqrt{\tau }L_1/L_2}\right){\displaystyle \frac{1}{2}}\mathrm{ln}L_2{\displaystyle \frac{1}{4}}\mathrm{ln}x_1{\displaystyle \frac{\pi \sqrt{\tau }L_1}{12L_2}}.`$
Introducing the Jacobi theta function

$$\vartheta _1(\varphi ,q)=2q^{1/4}\mathrm{sin}\varphi \underset{n=1}{\overset{\mathrm{}}{}}(1q^{2n})(12q^{2n}\mathrm{cos}2\varphi +q^{4n}),$$
(42)
and the identity

$$\underset{n=1}{\overset{\mathrm{}}{}}\left(1q^{2n}\right)=\left[\frac{\vartheta _1^{}(0,q)}{2q^{1/4}}\right]^{1/3}$$
(43)
where $`\vartheta _1^{}`$ is the derivative of $`\vartheta _1`$ with respect to $`\varphi `$, then we have also

$$c_3(x_1,x_2)=\frac{1}{2}\mathrm{ln}L_2+\frac{1}{4}\mathrm{ln}\frac{(x_1+x_2)}{x_1^3x_2^2}+\frac{1}{3}\mathrm{ln}\left[4\vartheta _1^{}(0,q)\right]$$
(44)
where

$$q=e^{L_1\sqrt{\tau }\pi /L_2}.$$
(45)
The identity (41) follows from the Jacobi transformation

$$\vartheta _1^{}(0,e^{\pi v})=v^{3/2}\vartheta _1^{}(0,e^{\pi /v}),v>0.$$
(46)
Now, the theta function $`\vartheta _1^{}`$ is finite for $`L_1/L_2`$ finite,<sup>1</sup><sup>1</sup>1This also follows from the fact that the series $`_{n=1}^{\mathrm{}}\mathrm{ln}(1u^n)=_{n=1}^{\mathrm{}}_{m=1}^{\mathrm{}}u^{mn}/m=_{m=1}^{\mathrm{}}u^m/m(1u^m)`$ converges for all $`0u<1`$. it follows that the leading behavior of $`c_3`$ is $`\mathrm{ln}L_2(\mathrm{ln}L_1)`$. In conformal field theory  one needs to compute the limits

$`{\displaystyle \frac{1}{L_1}}\underset{L_2\mathrm{}}{lim}{\displaystyle \frac{\mathrm{ln}T(𝐙_2)}{L_2}}`$ $`=`$ $`f_{\mathrm{bulk}}+{\displaystyle \frac{c_1}{L_1}}+{\displaystyle \frac{\mathrm{\Delta }_1}{L_1^2}}+o(L_1^2);`$ (47)
$`{\displaystyle \frac{1}{L_2}}\underset{L_1\mathrm{}}{lim}{\displaystyle \frac{\mathrm{ln}T(𝐙_2)}{L_1}}`$ $`=`$ $`f_{\mathrm{bulk}}+{\displaystyle \frac{c_2}{L_2}}+{\displaystyle \frac{\mathrm{\Delta }_2}{L_2^2}}+o(L_2^2).`$ (48)
Using (39), we find

$$\mathrm{\Delta }_1=\frac{\pi }{12\sqrt{\tau }},\mathrm{\Delta }_2=\frac{\pi \sqrt{\tau }}{12}.$$
(49)
This yields a central charge $`c=2`$ upon taking $`x_1=x_2=1`$ ($`\tau =1`$).

### 4.2 Dimer enumerations

We are now in a position to analyze the finite-size corrections of the two dimer enumerations (14) and (15). Although the expressions refer to two dimer lattices with different geometry, one for an $`M\times N1`$ lattice with a vacancy and $`MN=`$ odd, and one for an $`M\times N`$ lattice with $`MN=`$ even, a comparison can still be meaningful if both expansions are expressed in terms of lattice sizes $`M`$ and $`N`$.

(a) Close-packed dimers: For close-packed dimers on an $`M\times N`$ net with $`MN=`$ even, we have carried out the analysis for the expression (15) along the lines outlined in the above, and obtained the result (which can also be extracted from discussions in )

$$\mathrm{ln}Z_{\{M\times N\}}(x_1,x_2)=(MN+1)\overline{f}_{\mathrm{bulk}}+N\overline{c}_1+M\overline{c}_2+\overline{c}_3+o(1),$$
(50)
where

$`\overline{f}_{\mathrm{bulk}}(x_1,x_2)`$ $`=`$ $`{\displaystyle \frac{1}{4}}f_{\mathrm{bulk}}(x_1^2,x_2^2),`$
$`\overline{c}_1(x_1,x_2)`$ $`=`$ $`\overline{f}_{\mathrm{bulk}}(x_1,x_2){\displaystyle \frac{1}{2}}\mathrm{ln}\left(x_1+\sqrt{x_1^2+x_2^2}\right)`$
$`\overline{c}_2(x_1,x_2)`$ $`=`$ $`\overline{f}_{\mathrm{bulk}}(x_1,x_2){\displaystyle \frac{1}{2}}\mathrm{ln}\left(x_2+\sqrt{x_1^2+x_2^2}\right)`$
$`\overline{c}_3(x_1,x_2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}2{\displaystyle \frac{1}{2}}\mathrm{ln}\left(x_1+\sqrt{x_1^2+x_2^2}\right){\displaystyle \frac{1}{2}}\mathrm{ln}\left(x_2+\sqrt{x_1^2+x_2^2}\right)`$
$`+{\displaystyle \frac{1}{4}}\mathrm{ln}(x_1^2+x_2^2)+{\displaystyle \frac{\pi Mx_2}{24Nx_1}}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{ln}\left(1+e^{(2n1)\pi Mx_2/Nx_1}\right).`$
Expressions of $`\overline{c}_1(x_1,x_2)`$ and $`\overline{c}_2(x_1,x_2)`$ reduce to those found by Kenyon when $`x_1=x_2=1`$. In the language of the conformal field theory, the term $`\pi Mx_2/24Nx_1`$ in $`\overline{c}_3`$ yields the central charge $`c=1`$ upon taken $`M=N`$ and $`x_1=x_2`$, the accepted value for dimer and Ising systems.

Again, the expression (LABEL:dimercorrection1) for $`\overline{c}_3`$ is symmetric in $`\{x_1,M\}\{x_2,N\}`$, a fact can be seen from the identity

$$\frac{\pi Mx_2}{24Nx_1}+\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{ln}\left(1+e^{(2n1)\pi Mx_2/Nx_1}\right)=\frac{\pi Nx_1}{24Mx_2}+\underset{m=1}{\overset{\mathrm{}}{}}\mathrm{ln}\left(1+e^{(2m1)\pi Nx_1/Mx_2}\right).$$
(52)
The series $`_{n=1}^{\mathrm{}}\mathrm{ln}(1+u^{2n1})`$ converges,<sup>2</sup><sup>2</sup>2This also follows from the fact that the series $`_{n=1}^{\mathrm{}}\mathrm{ln}(1+u^{2n1})=_{m=1}^{\mathrm{}}(u)^m/m(1u^{2m})`$ converges for all $`0u<1`$. so $`\overline{c}_3`$ does not diverge for large $`M,N`$.

The expression for $`\overline{c}_3`$ can also be written as

$`\overline{c}_3(x_1,x_2)`$ $`=`$ $`{\displaystyle \frac{2}{3}}\mathrm{ln}2{\displaystyle \frac{1}{2}}\mathrm{ln}(x_1+\sqrt{x_1^2+x_2^2}){\displaystyle \frac{1}{2}}\mathrm{ln}(x_2+\sqrt{x_1^2+x_2^2})`$ (53)
$`+{\displaystyle \frac{1}{4}}\mathrm{ln}(x_1^2+x_2^2)+{\displaystyle \frac{1}{2}}\mathrm{ln}\vartheta _3(0,q){\displaystyle \frac{1}{6}}\mathrm{ln}\vartheta _1^{}(0,q),`$
where $`q=e^{\pi Mx_2/Nx_1}`$ and the theta function $`\vartheta _3`$ is given by

$$\vartheta _3(\varphi ,q)=\underset{n=1}{\overset{\mathrm{}}{}}(1q^{2n})(1+2q^{2n1}\mathrm{cos}2\varphi +q^{4n2}).$$
(54)
Then, the identity (52) is a consequence of the Jacobi transformations (41) and

$$\vartheta _3(0,e^{\pi /v})=v^{1/2}\vartheta _3(0,e^{\pi v}),v>0.$$
(55)
(b) Close-packed dimers with a boundary vacancy: For the simple-quartic net $`M\times N1`$ with a boundary vacancy and $`MN=`$ odd, one uses

$$L_1=(M+1)/2,L_2=(N+1)/2,$$
(56)
and expand (39) for large $`M`$ and $`N`$. After some algebra, we find

$`\mathrm{ln}Z_{\{M\times N1\}}(x_1,x_2)`$ (57)
$`=`$ $`(L_1L_2L_1)\mathrm{ln}x_1+(L_1L_2L_2)\mathrm{ln}x_2+\mathrm{ln}T(G;{\displaystyle \frac{x_1}{x_2}},{\displaystyle \frac{x_2}{x_1}})`$
$`=`$ $`(L_1L_2L_1)\mathrm{ln}x_1+(L_1L_2L_2)\mathrm{ln}x_2+L_1L_2\left[f_{\mathrm{bulk}}({\displaystyle \frac{x_1}{x_2}},{\displaystyle \frac{x_2}{x_1}})+f_c({\displaystyle \frac{x_1}{x_2}},{\displaystyle \frac{x_2}{x_1}})\right]`$
$`=`$ $`(MN+1)\overline{f}_{\mathrm{bulk}}+N\overline{c}_1+M\overline{c}_2+\overline{c}_3^{}+o(1),`$
where $`\overline{c}_1`$ and $`\overline{c}_2`$ are given in (LABEL:dimercorrection1), and

$`\overline{c}_3^{}(x_1,x_2)`$ $`=`$ $`{\displaystyle \frac{3}{2}}\mathrm{ln}{\displaystyle \frac{1}{2}}\mathrm{ln}\left(x_1+\sqrt{x_1^2+x_2^2}\right){\displaystyle \frac{1}{2}}\mathrm{ln}\left(x_2+\sqrt{x_1^2+x_2^2}\right)+{\displaystyle \frac{1}{4}}\mathrm{ln}\left(1+{\displaystyle \frac{x_2^2}{x_1^2}}\right)`$ (58)
$`{\displaystyle \frac{1}{2}}\mathrm{ln}N{\displaystyle \frac{\pi Mx_2}{12Nx_1}}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{ln}\left(1e^{2n\pi Mx_2/Nx_1}\right).`$
Comparing $`\overline{c}_3^{}`$ with $`\overline{c}_3`$ given by (LABEL:dimercorrection1), we see that for large $`M,N`$ the deletion of a boundary site introduces a logarithmic correction term $`\mathrm{ln}\sqrt{N}`$ in $`\overline{c}_3^{}`$. Furthermore, upon taking $`M=N`$ the term $`\pi Mx_2/12Nx_1`$ in $`\overline{c}_3^{}`$ yields a central charge $`c=2`$, which is different from that of the dimer system without vacancies.

To verify the occurrence of a logarithmic term for the vacancy problem, we consider the following example. Consider two dimer nets of $`N^21`$ sites each, where $`N3`$ is an odd integer so that the nets admit dimer coverings. While the two nets have different geometries, one a rectangular net of size $`(N+1)\times (N1)`$ and one a square net of size $`N\times N`$ with one boundary odd site removed, they have the same area and perimeter. Any difference in the evaluationsof (50) and (57) would occur in $`\overline{c}_3`$ and $`\overline{c}_3^{}`$ and higher order terms. Now from (50) and (57) we obtain

$`\mathrm{ln}Z_{\{(N+1)\times (N1)\}}(1,1)`$ $`=`$ $`N^2\left({\displaystyle \frac{G}{\pi }}\right)+2N\overline{c}_1\mathrm{ln}(1+\sqrt{2})+{\displaystyle \frac{3}{4}}\mathrm{ln}2+{\displaystyle \frac{\pi }{24}}`$
$`+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{ln}\left[1+e^{(2n1)\pi }\right]+o(1)`$
$`\mathrm{ln}Z_{\{N\times N1\}}(1,1)`$ $`=`$ $`(N^2+1)\left({\displaystyle \frac{G}{\pi }}\right)+2N\overline{c}_1\mathrm{ln}(1+\sqrt{2})+{\displaystyle \frac{7}{4}}\mathrm{ln}2`$ (59)
$`{\displaystyle \frac{\pi }{12}}{\displaystyle \frac{1}{2}}\mathrm{ln}N+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{ln}\left[1e^{2n\pi }\right]+o(1).`$
Defining the ratio

$$R(N)\frac{Z_{\{(N+1)\times (N1)\}}(1,1)}{Z_{\{N\times N1\}}(1,1)}$$
(60)
and using (59), we find the large $`N`$ behavior

$$R(N)C\sqrt{N},N\mathrm{},$$
(61)
where

$`C=\underset{N\mathrm{}}{lim}{\displaystyle \frac{R(N)}{\sqrt{N}}}`$ $`=`$ $`{\displaystyle \frac{e^{\pi /8G/\pi }}{2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1+e^{(2n1)\pi }}{1e^{2n\pi }}}\right).`$ (62)
$`=`$ $`\mathrm{0.\; 578\; 250}\mathrm{}.`$
As a numerical check, we have computed the value of $`R(N)/\sqrt{N}`$ for $`N=3`$ to $`2251`$ using (14) and (15). For $`N=9`$, for example, one has

$`Z_{\{10\times 8\}}(1,1)`$ $`=`$ $`\mathrm{1\; 031\; 151\; 241}=(89)\times (\mathrm{11\; 585\; 969})`$
$`Z_{\{9\times 91\}}(1,1)`$ $`=`$ $`\mathrm{557\; 568\; 000}=2^{12}\times 3^2\times 5^3\times (11)^2`$
$`R(9)/\sqrt{9}`$ $`=`$ $`\mathrm{0.616\; 457}\mathrm{}.`$ (63)
Results plotted in Fig. 3 confirm the large $`N`$ limit of $`C`$ given by (62) as well as the occurrence of the logarithmic correction in the vacancy problem. We remark that a similar result obtained by Kenyon  involving the occurrence of vacancy sites in the middle of a rectilinear net gives the ratio

$$R(N)N^{3/4},N\mathrm{},$$
(64)
and thus a logarithmic correction $`\frac{3}{4}\mathrm{ln}N`$ in the finite-size expansion.

## 5 Summary and discussions

We have used the Temperley-Kenyon-Propp-Wilson bijection between spanning trees on a lattice with free boundaries and dimer configurations on a related lattice with a boundary vacancy to establish the independence of the dimer generating function on the location of the vacancy. The equivalence is stated in Proposition (3). The generating function for close-packed dimers on a lattice with a single boundary vacancy is next computed, and compared with that of the known results for dimers without vacancies. It is found that the vacancy introduces a logarithmic correction in the large lattice expansion. A concrete example exhibiting this correction for an $`MN=`$ odd net with a vacancy as compared to an $`MN=`$ even net without vacancies is given.

To ascertain whether the logarithmic correction is due to the defect of a vacancy, or due to the oddness of the net size, one needs to compare expansions for two nets (of the same even-even lattice), one with two boundary vacancies and one without vacancies. While this problem can be formulated as the evaluation of the inverse of a matrix  in the Kasteleyn formulation, we argue that since the correction in question is that of the physical free energy of a dimer system, on physical ground one expects the correction to be additive for vacancies located sufficiently far apart. This would imply the logarithmic correction to be a “local” property due to the occurrence of vacancies.

We have also found that in the language of the conformal field theory the central charge for the vacancy problem is $`c=2`$ as compared to the value of $`c=1`$ for the dimer solution without vacancies. Furthermore, the $`\sqrt{N}`$ ratio (61) implies the existence of a boundary operator with scaling dimension 1/2, a value which does not appear in the standard Kac classification of operators at central charge $`2`$.<sup>3</sup><sup>3</sup>3We are indebted to the referee for this remark. The extraction of the central charge should be viewed with caution, however, since the dimer systems do not exhibit critical points.

We have also carried out finite-size analyses (details of which to be given elsewhere) for spanning trees on simple-quartic nets with other, including the toroidal, cylindrical, Möbius, and Klein bottle, boundary conditions. It is found that the logarithmic correction reported in this paper arises only in the case of free boundaries. This is consistent to the fact that the formulation of the extended Temperley bijection as presented in this paper is a property that is unique to graphs with free boundaries.

## Acknowledgment

Work has been supported in part by the National Science Council of the Republic of China under Grant No. NSC 90-2112-M-032-003 and the National Science Foundation Grant DMR-9980440. We thank W. T. Lu for a useful discussion and R. Kenyon for calling our attention to . We also thank the referee for suggesting to us the formulation of Proposition 3.1.2.

Figure captions

Fig. 1. The construction of dimer lattices $`G_D`$ from a $`3\times 3`$ spanning tree lattice $`G`$. Solid circles denote odd sites and open circles denote the odd site that has been removed in $`G_D`$. (a) A spanning tree lattice $`G`$. (b) A dimer lattice $`G_D`$ constructed from $`G`$ with one corner site removed. (c) A dimer lattice $`G_D`$ constructed from $`G`$ with one odd site on the boundary removed.

Fig. 2. The bijection between spanning trees on $`G`$ and dimer configurations of $`G_D`$. (a) A spanning tree configuration on $`G`$. (b) The corresponding dimer configuration on the dimer lattice $`G_D`$ of Fig. 1(b). (c) The corresponding dimer configuration on the dimer lattice $`G_D`$ of Fig. 1(c).

Fig. 3. The enumeration of $`R(N)/\sqrt{N}`$ for $`N=`$ 3 to 2251. The dashed line indicates the value $`C`$ given by (62) in the large $`N`$ limit.


format: text

text:
# 1 Introduction

## 1 Introduction

Experiments at hadron colliders require apart from the hard scattering amplitudes, which can be calculated in perturbation theory, also the knowledge of the parton distributions functions (pdf’s). Although these parton distribution functions parameterize the non-perturbative information, they depend on a factorization scale and the dependence on that scale can be calculated within perturbation theory. Usually parameterizations of parton distributions functions are given at a low input scale  \-  and then evolved to the desired scale. At present the evolution kernel is known completely to next-to-leading order (NLO)  \- . The evolution has to be done numerically. There are two methods available: The first one consists in integrating numerically the evolution equations in $`x`$ and in $`Q^2`$ and is called the x-space method. The second one first performs a Mellin-transform on the variable $`x`$ . The evolution equations then factorize and can be solved analytically. Using this Mellin-transform method (or N-space method), leaves us with the task to evaluate the inverse Mellin transform numerically. In this paper I describe a numerical program, which uses an optimized integration contour in the complex N-plane in order to do the inverse Mellin transform. This method has been proposed by Kosower . Using this method gives a fast and accurate program for the evolution of parton distributions. Although a variety of numerical evolution programs already exist (mainly within the x-space method)  \- , they usually cannot compete in computation speed. The required CPU time becomes an issue if one tries to extend the analysis of experiments at hadron colliders in two important directions:

* Parton distributions with error-bars . The Tevatron data on the measurement of the one jet inclusive transverse energy has shown  that the method of obtaining parton distributions from “global fits” has reached its limits. Giele, Keller and Kosower  have proposed a method to include uncertainties into the parton distributions. This method involves a functional integration over a set of parton distributions. With a set of pdf’s containing 1000 or more parton distributions, the traditional approach of a pre-calculated grid containing the evolved parton distributions becomes inappropriate. Instead, one would need a program which allows to calculate the evolution on the fly. The program presented here allows one to do just that.
* Extension to NNLO. The calculation of the three-loop anomalous dimensions is currently under way . It is likely that the result will be rather lengthy. The computational cost of evaluating the anomalous dimensions will therefore be rather high, and one would like to minimize the number of function calls to the anomalous dimensions. The optimized contour employed in this program minimizes this number. Once the three-loop anomalous dimensions are known, the program can be extended in a straightforward way to include NNLO corrections.

The program is optimized for the N-space method, where one needs typically only four or five evaluations of the anomalous dimensions to obtain the evolution of a parton distribution with an accuracy of $`210^4`$. In addition, I also include an option to simulate the x-space evolution programs. x-space programs use a different truncation prescription and the results might differ numerically from the N-space approach. It should be noted, that both approaches are valid to next-to-leading order, the difference is formally of higher order. The x-space option is included to allow an easy comparison with those programs, however it is not optimized. In technical terms, since no closed formula for the evolution operator in the singlet sector within the x-space approach is known, the evolution operator is approximated by a sequence of evolutions over a small interval. This slows down the program considerably. The N-space evolution method is therefore the recommended one, where the full performance of the program is obtained.
This paper is organized as follows. In section 2 I recall the basic formulae for the evolution of parton densities and sketch how to choose an optimized contour using the Mellin transform method. Section 3 gives an overview of the design of the program, while section 4 is of a more practical nature and describes how to install and use the program. In section 5 I compare the program with the reference results of Blümlein et al.  and address issues like performance. A summary is provided in section 6. The appendix collects the formulae for the evolution operators.

## 2 Theoretical Background

The evolution equations for the parton distributions $`f(x,Q^2)`$ of the proton read

$`Q^2{\displaystyle \frac{f(x,Q^2)}{Q^2}}`$ $`=`$ $`P(x,Q^2)f(x,Q^2),`$ (1)
where $`x`$ stands for the nucleon’s momentum fraction carried by the parton, $`P(x,Q^2)`$ is the Altarelli-Parisi evolution kernel, and $``$ denotes the convolution

$`A(x)B(x)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{1}{}}}𝑑y{\displaystyle \underset{0}{\overset{1}{}}}𝑑z\delta (xyz)A(y)B(z).`$ (2)
Eq. (1) is understood to represent the quark non-singlet evolution as well as the coupled quark singlet and gluon evolution. In the latter case $`f`$ is a two-component vector and $`P(x,Q^2)`$ a two-by-two matrix. Eq. (1) can be factorized by taking Mellin moments

$`A^z`$ $`=`$ $`{\displaystyle \underset{0}{\overset{1}{}}}𝑑xx^{z1}A(x)`$ (3)
and one obtains

$`Q^2{\displaystyle \frac{f^z(Q^2)}{Q^2}}`$ $`=`$ $`P^z(Q^2)f^z(Q^2).`$ (4)
Changing variables by using $`a_s=\alpha _s(Q^2)/4\pi `$ instead of $`\mathrm{ln}Q^2`$ as evolution variable one obtains

$`{\displaystyle \frac{f^z(Q^2)}{a_s}}`$ $`=`$ $`{\displaystyle \frac{P^z(a_s)}{\beta (a_s)}}f^z(Q^2),`$ (5)
where

$`\beta (a_s)`$ $`=`$ $`Q^2{\displaystyle \frac{a_s}{Q^2}}`$ (6)
is the beta function for the strong coupling. Up to now, the evolution kernel is known completely only to next-to-leading order:

$`P^z(Q^2)`$ $`=`$ $`a_sP_0^z+a_s^2P_1^z+O(a_s^3).`$ (7)
To this order, the beta function is given by

$`\beta (a_s)`$ $`=`$ $`\beta _0a_s^2\beta _1a_s^3+O(a_s^4).`$ (8)
Inserting eq. (7) and eq. (8) into eq. (5) and expanding the r.h.s in $`a_s`$ one obtains

$`{\displaystyle \frac{f^z(Q^2)}{a_s}}`$ $`=`$ $`{\displaystyle \frac{1}{\beta _0a_s}}\left[P_0^z+a_s\left(P_1^z{\displaystyle \frac{\beta _1}{\beta _0}}P_0^z\right)\right]f^z(Q^2).`$ (9)
Solutions of eq. (9) are referred to as “N-space” solutions, since this is the standard approach within the N-space method. One the other hand, one might defer from inverting the power series in the denominator on the r.h.s of eq. (5) and to solve eq. (5) directly. To NLO, this amounts to solving

$`{\displaystyle \frac{f^z(Q^2)}{a_s}}`$ $`=`$ $`{\displaystyle \frac{\left(P_0^z+a_sP_1^z\right)}{a_s\left(\beta _0+\beta _1a_s\right)}}f^z(Q).`$ (10)
Solutions of eq. (10) are referred to as “x-space” solutions, again since this corresponds to the usual procedure within the x-space method. Eq. (9) and eq. (10) differ by terms which are formally of higher order in $`a_s`$, and solutions to both equations are therefore considered accurate to next-to-leading order, although they might differ numerically.
In addition there are different common practices on how to calculate $`a_s`$ at the scale $`Q^2`$. At NLO $`a_s(Q^2)`$ can be found from the solution of the equation

$`{\displaystyle \frac{1}{a_s}}{\displaystyle \frac{\beta _1}{\beta _0}}\mathrm{ln}\left(\beta _1+{\displaystyle \frac{\beta _0}{a_s}}\right)\beta _0L`$ $`=`$ $`0,`$ (11)
where $`L=\mathrm{ln}\frac{Q^2}{\mathrm{\Lambda }_{QCD}^2}`$. This equation has to be solved numerically, and solutions of eq. (11) are referred to as “exact” NLO solutions. On the other hand it is common practice to view $`a_s`$ as a power series in $`1/L`$ and to use

$`a_s`$ $`=`$ $`{\displaystyle \frac{1}{\beta _0L}}\left(1{\displaystyle \frac{\beta _1}{\beta _0^2}}{\displaystyle \frac{\mathrm{ln}L}{L}}\right)`$ (12)
to calculate $`a_s`$ at the scale $`Q^2`$. Solutions of eq. (12) are referred to as “truncated in $`1/L`$” solutions. Again, the numerical values might differ for solutions of eq. (11) and eq. (12). It should be mentioned that the integration constant $`\mathrm{\Lambda }_{QCD}`$ appearing in eq. (11) is not identical to the integration constant of the same name in eq. (12), e.g. fixing $`a_s`$ at a scale like $`m_Z^2`$ and solving eq. (11) and eq. (12) for $`\mathrm{\Lambda }_{QCD}`$ will yield two different numerical values. Each value has to be used in conjunction with the equation from which it was obtained.
Eq. (9) and eq. (10) are formally solved by the introduction of an evolution operator $`E^z`$:

$`f^z(Q^2)`$ $`=`$ $`E^z(a_s(Q^2),a_s(Q_0^2))f^z(Q_0^2)`$ (13)
Furthermore, there is yet another truncation method used in the literature: Within this method one does not perform the change of the evolution variable from $`Q^2`$ to $`a_s`$, but one substitutes the truncated solution eq. (12) for $`a_s`$ into eq. (4). This amounts to solving

$`Q^2{\displaystyle \frac{f^z(Q^2)}{Q^2}}`$ $`=`$ $`\left({\displaystyle \frac{1}{\beta _0L}}\left(1{\displaystyle \frac{\beta _1}{\beta _0^2}}{\displaystyle \frac{\mathrm{ln}L}{L}}\right)P_0^z+{\displaystyle \frac{1}{\beta _0^2L^2}}\left(1{\displaystyle \frac{\beta _1}{\beta _0^2}}{\displaystyle \frac{\mathrm{ln}L}{L}}\right)^2P_1^z\right)f^z(Q^2).`$ (14)
I will refer to solutions of eq. (14) as “x-space and truncated in $`1/L`$” solutions. Again, this equation is formally solved by an evolution operator, which now depends on $`Q^2`$ and $`Q_0^2`$ instead of $`a_s(Q^2)`$ and $`a_s(Q_0^2)`$:

$`f^z(Q^2)`$ $`=`$ $`\widehat{E}^z(Q^2,Q_0^2)f^z(Q_0^2)`$ (15)
Formulae for the evolution operators corresponding to the various truncation prescriptions are collected in the appendix. The evolved parton distributions in x-space are then obtained by an inverse Mellin transformation:

$`f(x,Q^2)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle \underset{C}{}}𝑑zx^zf^z(Q^2).`$ (16)
Here, the contour $`C`$ runs to the right of all singularities of the integrand.
The parton distributions are usually parametrized at the input scale $`Q_0^2`$ in a form like

$`xf(x,Q_0^2)`$ $`=`$ $`{\displaystyle \underset{i}{}}A_ix^{\alpha _i}(1x)^{\beta _i}`$ (17)
with Mellin transform

$`f^z(Q_0^2)`$ $`=`$ $`{\displaystyle \underset{i}{}}A_iB(z+\alpha _i1,1+\beta _i),`$ (18)
where $`B(x,y)`$ is Euler’s beta function. Dropping from now on the arguments $`Q^2`$ and $`Q_0^2`$ our task is to evaluate the integral

$`I`$ $`=`$ $`\text{Re}{\displaystyle \frac{1}{\pi i}}{\displaystyle \underset{C_s}{}}𝑑zE^zF(z),`$ (19)
$`F(z)=x^z{\displaystyle \underset{i}{}}A_iB(z+\alpha _i1,1+\beta _i),`$
where complex conjugation has been used to replace the contour $`C`$ by $`C_s`$, starting at the real axis, right to the right-most pole and running upwards to infinity. In the singlet case $`F(z)`$ and $`A_i`$ are vectors and $`E^z`$ is a matrix. The most efficient way to solve this problem is to choose a contour in such a way that the integrand can very well be approximated by some set of orthogonal polynomials. The method is due to Kosower . It uses a parabolic contour, completely determined by the parton distribution at the original scale $`Q_0^2`$ and evaluates the integral by a Gauss-Laguerre quadrature formula. With this method a few (e.g. four or five) evaluations of the integrand are sufficient to evaluate the integral to very high precision (e.g. to an accuracy of 0.02 %).
I shortly review this method for the non-singlet case. The contour of integration is given by

$`z(u)=z_0+ic_2\sqrt{u}+{\displaystyle \frac{1}{2}}c_2^2c_3u,u=0\mathrm{}\mathrm{},`$ (20)
where $`z_0`$ is the minima of the function $`F(z)`$ on the real axis right to the rightmost pole of $`F(z)`$. The parameters $`c_2`$ and $`c_3`$ are given by

$`c_2=\sqrt{{\displaystyle \frac{2F(z_0)}{F^{\prime \prime }(z_0)}}},`$ $`c_3={\displaystyle \frac{F^{\prime \prime \prime }(z_0)}{3F^{\prime \prime }(z_0)}}.`$ (21)
$`F^{}(z)`$, $`F^{\prime \prime }(z)`$ and $`F^{\prime \prime \prime }(z)`$ are the first three derivatives of $`F(z)`$ with respect to $`z`$. For the parameterization eq. (17) they are given by

$`F^{}`$ $`=`$ $`{\displaystyle \underset{i}{}}G_i^{}F_i,`$
$`F^{\prime \prime }`$ $`=`$ $`{\displaystyle \underset{i}{}}\left(G_i^{\prime \prime }+G_i^2\right)F_i,`$
$`F^{\prime \prime \prime }`$ $`=`$ $`{\displaystyle \underset{i}{}}\left(G_i^{\prime \prime \prime }+3G_i^{\prime \prime }G_i^{}+G_i^3\right)F_i,`$ (22)
where $`F_i=x^zA_iB(z+\alpha _i1,1+\beta _i)`$ and

$`G_i`$ $`=`$ $`\mathrm{ln}(x^zB(z+\alpha _i1,1+\beta _i)),`$
$`G_i^{}`$ $`=`$ $`\mathrm{ln}x+\psi (z+\alpha _i1)\psi (z+\alpha _i+\beta _i),`$
$`G_i^{\prime \prime }`$ $`=`$ $`\psi ^{}(z+\alpha _i1)\psi ^{}(z+\alpha _i+\beta _i),`$
$`G_i^{\prime \prime \prime }`$ $`=`$ $`\psi ^{\prime \prime }(z+\alpha _i1)\psi ^{\prime \prime }(z+\alpha _i+\beta _i).`$ (23)
Here, $`\psi `$, $`\psi ^{}`$ and $`\psi ^{\prime \prime }`$ are polygamma functions. With this parameterization one obtains for the integral eq. (19):

$`I`$ $`=`$ $`{\displaystyle \frac{c_2}{2\pi }}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{du}{\sqrt{u}}}e^u\text{Re}\left[e^u\left(1ic_2c_3\sqrt{u}\right)E^{z(u)}F\left(z(u)\right)\right].`$ (24)
The function $`u^{1/2}e^u`$ is the weight function of the Laguerre polynomials $`L_k^{1/2}(x)`$ and the integral can be approximated by a Gauss-Laguerre quadrature formula:

$`I`$ $``$ $`{\displaystyle \frac{c_2}{2\pi }}{\displaystyle \underset{j=1}{\overset{k}{}}}w_j\text{Re}\left[e^{u_j}\left(1ic_2c_3\sqrt{u_j}\right)E^{z(u_j)}F\left(z(u_j)\right)\right].`$ (25)
Here, $`u_j`$ are the zeros of $`L_k^{1/2}(x)`$ and the weights $`w_j`$ are given by

$`w_j`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }\left(n+\frac{1}{2}\right)}{n!(n+1)^2}}{\displaystyle \frac{u_j}{\left(L_{k+1}^{1/2}(u_j)\right)^2}}.`$ (26)
In the singlet case one has a two-component vector given by the quark singlet combination $`\mathrm{\Sigma }^z`$ and the gluon distribution $`G^z`$. There will be strong mixing already for moderate evolution in $`Q^2`$, since $`\mathrm{\Sigma }^z`$ and $`G^z`$ are not eigenfunctions of the evolution kernel. One might be tempted to work in a basis of eigenfunctions of the leading-order evolution kernel and to extend the method discussed above to these eigenfunctions, eventually choosing two different contours for the two eigenfunctions. However this approach faces two problems: First of all, one of the two eigenfunctions usually does not have a minimum on the real axis on the right of the rightmost pole. Secondly, by using different contours for the two components one is forced to keep track of branch cut crossings, a task one would like to avoid. It is therefore better to use a single contour for both components, although it might not be the optimal one for both components.
In the singlet sector the contour is chosen as follows (this differs from the original proposition of Kosower): $`z_0`$ is chosen as the minimum of the sum of the $`\mathrm{\Sigma }^z`$\- and the $`G^z`$-component. I then diagonalize the leading-order evolution kernel at $`z_0`$:

$`\left(\begin{array}{c}e_1^z\\ e_2^z\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\gamma _{0,qg}(z_0)& \gamma _{0,qq}(z_0)\lambda _{}(z_0)\\ \gamma _{0,gg}(z_0)\lambda _+(z_0)& \gamma _{0,gq}(z_0)\end{array}\right)\left(\begin{array}{c}\mathrm{\Sigma }^z\\ G^z\end{array}\right),`$ (33)
where the eigenvalues $`\lambda _\pm `$ are given by

$`\lambda _\pm (z_0)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\gamma _{qq}^{(0)}(z_0)+\gamma _{gg}^{(0)}(z_0)\pm \sqrt{\left(\gamma _{gg}^{(0)}(z_0)\gamma _{qq}^{(0)}(z_0)\right)^2+4\gamma _{qg}^{(0)}(z_0)\gamma _{gq}^{(0)}(z_0)}\right].`$ (34)
Note that any $`z`$-dependence in eq. (33) is only due to $`\mathrm{\Sigma }^z`$ or $`G^z`$, $`z_0`$ appearing in the anomalous dimensions is a fixed parameter. Strictly speaking the eigenvectors of eq. (33) are only eigenvectors at the point $`z=z_0`$ of the leading-order evolution kernel. I then take $`e_1^z`$ to determine $`c_2`$ and $`c_3`$ and integrate both components with this contour. To compensate for the fact that $`z_0`$ is not necessarily the minimum of $`e_1^z`$, the formula eq. (21) is replaced by

$`c_2`$ $`=`$ $`\sqrt{{\displaystyle \frac{2F(z_0)F^{\prime \prime }(z_0)}{F^{\prime \prime }(z_0)^2F^{}(z_0)F^{\prime \prime \prime }(z_0)}}},`$ (35)
which takes into account corrections due to the fact, that the contour does not start from a minimum of $`F(z)`$ (e.g. $`F^{}(z_0)0`$). The reasons why this works is that in the cases of interest $`e_1^z`$ is roughly the sum of $`\mathrm{\Sigma }^z`$ and $`G^z`$ and so $`z_0`$ is approximately optimal for $`e_1^z`$. The other component $`e_2^z`$ is roughly the difference and numerically therefore smaller.

## 3 Design of the program

The program is written in C++. All classes and functions are defined in the namespace `pdf`. The parameterization of a parton distribution at the initial scale $`Q_0^2`$ is represented by the class `partondistribution`. A `partondistribution` is constructed as follows:

```
  // u-valence CTEQ 4M
  int n = 2;
  double Q0 = 1.6;
  double A_u[2]     = { 1.344, 1.344*6.402};
  double alpha_u[2] = { 0.501, 0.501+0.873 };
  double beta_u[2]  = { 3.689, 3.689 };
  int eta = -1;

  partondistribution uvalence =
    partondistribution(n,Q0,A_u,alpha_u,beta_u,eta);

```

This corresponds to the parameterization

$`xu_v(x,Q_0^2)`$ $`=`$ $`1.344x^{0.501}(1x)^{3.689}\left(1+6.402x^{0.873}\right)`$ (36)
at the initial scale $`Q_0=1.6\text{GeV}`$. Valence like distributions need to have $`\eta =1`$ explicitly in the constructor. For distributions which correspond to $`\eta =1`$ (e.g. distributions in the singlet sector or non-singlet non-valence distributions like $`(u+\overline{u})(d+\overline{d})`$), the $`\eta `$-parameter can be dropped.
This partondistribution can be evolved to a scale $`Q`$ by calling the function `evolve_nonsinglet`:

```
  // evolution to 10 GeV for x = 0.01
  double Q = 10.0;
  double x = 0.01;

  double u_v = evolve_nonsinglet(uvalence,Q,x);

```

The function `evolve_nonsinglet` is just an entry point and the job is passed to the function `evolve_nonsinglet_5`, which performs a Gauss-Laguerre quadrature with 5 points. In addition there are functions `evolve_nonsinglet_3`, `evolve_nonsinglet_10`,
`evolve_nonsinglet_20` and `evolve_nonsinglet_30`, corresponding to quadrature formulae with 3, 10, 20 and 30 points. These functions can be selected by passing an additional argument to `evolve_nonsinglet`, e.g.

```
  // evolution using a 10-point quadrature formula
  double u_v = evolve_nonsinglet(uvalence,Q,x,10);

```

will call the function `evolve_nonsinglet_10`.
In the singlet sector evolution is done with the help of the function `evolve_singlet`:

```
  // singlet evolution
  vector_d singlet = evolve_singlet(Sigma,Gluon,Q,x);

  double s = singlet[0];
  double g = singlet[1];

```

`Sigma` and `Gluon` are of type `partondistribution` and represent the parameterization of the quark singlet combination and the gluon distribution, respectively. The evolution function returns a two-component vector, with the evolved quark singlet distribution as first component and the evolved gluon distribution as second component. In C and C++ subscription of arrays starts at zero. As in the non-singlet case the function `evolve_singlet` is just an entry point, calling `evolve_singlet_5`. Again, there are functions using quadrature formulae with 3, 5, 10, 20 and 30 points.
In addition there are the classes `evolutionkernel` and `alpha_strong`. The class
`evolutionkernel` has a static data member `method`, which allows one to select the truncation method, e.g.

```
  // select truncation method for evolution kernel
  evolutionkernel::method = evolutionkernel::N_space;

```

will perform the evolution according to eq. (9). Other choices are `evolutionkernel::x_space` and `evolutionkernel::x_space_truncate_in_one_over_L`, corresponding to evolution according to eq. (10) and eq. (14), respectively.
In a similar way, the class `alpha_strong` has a static data member `method`, which allows one to select the truncation method for $`\alpha _s`$, e.g.

```
  // select truncation method for alpha_strong
  alpha_strong::method = alpha_strong::truncate_in_one_over_L;

```

will use formula eq.(12), whereas the choice `alpha_strong::exact` will use eq. (11). In addition the user should set the flavor thresholds and the values of $`\mathrm{\Lambda }_{QCD}`$ for the different number of flavors, as for example in

```
  // define values for alpha_strong (in GeV)
  alpha_strong::charm_threshold  = 1.3;
  alpha_strong::bottom_threshold = 4.7;

  alpha_strong::lambda_3 = 0.374;
  alpha_strong::lambda_4 = 0.327;
  alpha_strong::lambda_5 = 0.226;

```

Other classes of general interest provided by the package are the class `complex_d` corresponding to complex numbers with double precision as well as the templates `vector_template` `<class T>` and `matrix_template<class T>`. The specializations `vector_d`, `vector_c`,
`matrix_d` and `matrix_c` define vectors and matrices with real or complex entries in double precision.

## 4 How to Use the Library

In this section I give indications how to install and use the program library. Compilation of the package will build a (shared) library called `libpartonevolution`. The user will then link his own programs against the library.

### 4.1 Installation

The program library can be obtained from
http://www.fis.unipr.it/~stefanw/partonevolution
After unpacking, the library for the evolution of parton distributions is build by issuing the commands

```
  ./configure
  make
  make install

```

There are various options which can be passed to the configure script, an overview can be obtained with ./configure --help.
After installation, the shell script `partonevolution-config` can be used to determine the compiler and linker command line options required to compile and link a program with the partonevolution library. For example, `partonevolution-config --cppflags` will give the path to the header files of the library, whereas `partonevolution-config --libs` prints out the flags necesarry to link a program against the library.

### 4.2 Writing programs using the library

Once the library is build and installed, it is ready to be used. Here is a small example program, which defines a parameterization at the scale $`Q_0=2\text{GeV}`$ and calls the evolution routine for $`Q=10\text{GeV}`$ and $`x=0.01`$ .

```
#include <iostream>
#include "partonevolution/partonevolution.h"

int main()
{
  using namespace pdf;

  // select truncation methods
  alpha_strong::method = alpha_strong::truncate_in_one_over_L;
  evolutionkernel::method = evolutionkernel::N_space;

  // define a toy model with 4 flavors
  alpha_strong::charm_threshold  = 0.0;
  alpha_strong::bottom_threshold = 1E10;
  alpha_strong::lambda_4 = 0.250;

  // define a parton distribution
  double Q0 = 2.0;
  double A_u[1]     = { 35/16.};
  double alpha_u[1] = { 0.5 };
  double beta_u[1]  = { 3.0 };
  int eta = -1;

  partondistribution f =
    partondistribution(1,Q0,A_u,alpha_u,beta_u,eta);

  // evolution
  double Q = 10.0;
  double x = 0.01;

  double res = evolve_nonsinglet(f,Q,x);

  // print out result
  cout << "x*f(x,Q^2) = " << x*res << endl;
}

```

After compilation and linking against the partonevolution library, one obtains an executable, which will print out

```
x*f(x,Q^2) = 0.247242

```

This program defines the $`u_v`$-distribution

$`xu_v(x,Q_0^2)`$ $`=`$ $`{\displaystyle \frac{35}{16}}x^{0.5}(1x)^{3.0}`$ (37)
of the toy model of ref. . The numerical difference of our value with the value 0.24723 quoted in ref.  is well below the estimated accuracy given in ref. .

### 4.3 Documentation

The complete documentation of the program is inserted as comment lines in the source code. The documentation can be extracted from the sources with the help of the documentation system “doxygen” . The program “doxygen” is freely available. Issuing in the top-level build directory for the partonevolution library the commands

```
  doxygen Doxyfile

```

will create a directory “reference” with the documentation in html and latex format.

## 5 Checks and performance

In this section I compare the results of our program with the reference results of Blümlein et al. . The authors of ref.  define a toy model, consisting of four active flavors with $`\mathrm{\Lambda }_{QCD}^{(4)}=250\text{MeV}`$. At the input scale $`Q_0^2=4\text{GeV}^2`$ they take the initial parton distributions as follows:

$`xu_v=A_ux^{0.5}(1x)^3,`$ $`xd_v=A_dx^{0.5}(1x)^4,`$
$`xS=A_Sx^{0.2}(1x)^7,`$ $`xg=A_gx^{0.2}(1x)^5,`$
$`xc=0,`$ $`x\overline{c}=0.`$ (38)
The sea $`S`$ is taken to be $`SU(3)_{flavor}`$-symmetric and to carry $`15\%`$ of the nucleon’s momentum at the input scale. This, together with the usual flavor- and momentum sum rules, fixes the constants $`A_u`$, $`A_d`$, $`A_S`$ and $`A_g`$. The sea is related to the quark singlet distribution

$`\mathrm{\Sigma }={\displaystyle \underset{q=u,d,\mathrm{}}{}}\left(q+\overline{q}\right)`$ (39)
by

$`S`$ $`=`$ $`\mathrm{\Sigma }u_vd_v.`$ (40)
They provide a table with numerical results for the evolved parton distributions at the scale $`Q^2=100\text{GeV}^2`$. They have estimated the numerical accuracy of their results to be $`0.02\%`$. I have compared the program with the contents of table 1 of . Using a 5-point quadrature formula the results of the program agree with the numbers presented there within the indicated estimated accuracy (0.02%), with some exceptions for the charm distribution and the singlet case for small values of $`x`$. Using a 10-point quadrature formula for those cases will lead to results within the $`0.02\%`$ uncertainty margin. I have compared both the x-space solution (upper part of the table 1 in ref. , corresponding to solutions of eq. (14) ) as well as the N-space solution (lower part of table 1 in ref. , corresponding to solutions of eq. (9) in combination with eq. (12) for the strong coupling).
Finally, I would like to give some indication on the required CPU time. The evolution of a non-singlet parton distribution with a 5-point Gauss-Laguerre quadrature formula takes about $`1.310^3s`$ on a standard PC (Athlon, 1.6 GHz). The singlet evolution requires the evaluation of four anomalous dimensions. In addition there is some overhead for vector and matrix arithmetic. The singlet evolution is therefore more expensive in terms of CPU time. The evolution of the singlet combination takes about $`610^3s`$, again with a 5-point Gauss-Laguerre quadrature formula on the same PC. These numbers apply to the N-space method. Simulating the x-space evolution method, the evolution of a non-singlet parton distribution with a 5-point Gauss-Laguerre quadrature formula takes about $`1.310^3s`$. In the singlet sector, the evolution kernel has to be approximated by evolutions over small intervals. Subdividing the evolution interval into $`1000`$ subintervals, the evolution of the singlet combination takes about $`27010^3s`$, again with a 5-point Gauss-Laguerre quadrature formula on the same PC.

## 6 Summary

In this paper I described the program library “partonevolution”. This library can be used to evolve parton distributions given by parameterizations at a scale $`Q_0`$ to the desired scale $`Q`$. The program uses the Mellin transform method with an optimized contour. Due to this optimized contour the program needs only a few evaluations of the anomalous dimensions and is therefore extremely fast. The library is therefore suited to be used in analysis involving parton distributions with errors, where a functional integration over a set of pdf’s is needed. Furthermore, once the three-loop anomalous dimensions are known, the extension of the algorithms to NNLO is straightforward. In addition, the program can also be used to simulate x-space evolution programs.

## Appendix A Evolution operators

In this appendix the formulae for the evolution operators are collected. I start with the N-space method. The evolution operator, which solves eq. (9) is given in the non-singlet case by

$`E_\eta ^z(a_s(Q^2),a_s(Q_0^2))`$ $`=`$ $`\left({\displaystyle \frac{a_s(Q^2)}{a_s(Q_0^2)}}\right)^{\frac{\gamma _0^z}{2\beta _0}}\left[1+{\displaystyle \frac{a_s(Q^2)a_s(Q_0^2)}{2\beta _0}}\left(\gamma _1^z(\eta ){\displaystyle \frac{\beta _1}{\beta _0}}\gamma _0^z\right)\right].`$ (41)
Here, $`\eta =1`$ corresponds to the combinations $`q\overline{q}`$, while $`\eta =+1`$ should be used for the non-singlet combinations like $`(u+\overline{u})(d+\overline{d})`$. I have expressed the evolution operator in terms of the anomalous dimensions $`\gamma _0^z`$ and $`\gamma _1^z`$, which are related to the coefficients of the evolution kernel by $`P_i^z=\frac{1}{2}\gamma _i^z`$. The explicit expressions for the anomalous dimensions can be found in . For the singlet case one finds :

$`E^z(a_s(Q^2),a_s(Q_0^2))`$ $`=`$ $`\left({\displaystyle \frac{a_s(Q^2)}{a_s(Q_0^2)}}\right)^{\frac{\lambda _{}^z}{2\beta _0}}[P_{}^z+{\displaystyle \frac{a_s(Q^2)a_s(Q_0^2)}{2\beta _0}}P_{}^z\gamma ^zP_{}^z`$ (42)
$`(a_s(Q_0^2)a_s(Q^2)\left({\displaystyle \frac{a_s(Q^2)}{a_s(Q_0^2)}}\right)^{\frac{\lambda _+^z\lambda _{}^z}{2\beta _0}}){\displaystyle \frac{P_{}^z\gamma ^zP_+^z}{2\beta _0+\lambda _+^z\lambda _{}^z}}]`$
$`+(+)`$
with

$`\gamma ^z`$ $`=`$ $`\gamma _1^z{\displaystyle \frac{\beta _1}{\beta _0}}\gamma _0^z,`$
$`\lambda _\pm ^z`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\gamma _{0,qq}^z+\gamma _{0,gg}^z\pm \sqrt{\left(\gamma _{0,gg}^z\gamma _{0,qq}^z\right)^2+4\gamma _{0,qg}^z\gamma _{0,gq}^z}\right),`$
$`P_\pm ^z`$ $`=`$ $`\pm {\displaystyle \frac{1}{\lambda _+^z\lambda _{}^z}}\left(\gamma _0^z\lambda _{}^z\right).`$ (43)
The evolution operator solving eq. (10) (e.g. “x-space method”) reads in the non-singlet sector

$`E_\eta ^z(a_s(Q^2),a_s(Q_0^2))`$ $`=`$ $`\left({\displaystyle \frac{a_s(Q^2)}{a_s(Q_0^2)}}\right)^{\frac{\gamma _0^z}{2\beta _0}}\left({\displaystyle \frac{\beta _0+\beta _1a_s(Q^2)}{\beta _0+\beta _1a_s(Q_0^2)}}\right)^{\frac{\gamma ^z(\eta )}{2\beta _1}}.`$ (44)
In the singlet case no closed form is known and the evolution operator has to be obtained as a sequence of evolutions over a small interval:

$`E^z(a_s(Q^2),a_s(Q_0^2))`$ $`=`$ $`E^z(a^{(n)},a^{(n1)})E^z(a^{(n1)},a^{(n2)})\mathrm{}E^z(a^{(1)},a^{(0)}),`$ (45)
where $`a^{(j)}=a_s(Q_0^2)+\left(a_s(Q^2)a_s(Q_0^2)\right)j/n`$. In a small interval the evolution operator for the singlet case can be approximated by

$`E^z(a_s(Q^2),a_s(Q_0^2))`$ (46)
$`=`$ $`\left({\displaystyle \frac{a_s(Q^2)}{a_s(Q_0^2)}}\right)^{\frac{\lambda _{}^z}{2\beta _0}}[P_{}^z+{\displaystyle \frac{1}{2\beta _1}}\mathrm{ln}\left({\displaystyle \frac{\beta _0+\beta _1a_s(Q^2)}{\beta _0+\beta _1a_s(Q_0^2)}}\right)P_{}^z\gamma ^zP_{}^z+{\displaystyle \frac{P_{}^z\gamma ^zP_+^z}{2\beta _0+\lambda _+^z\lambda _{}^z}}`$
$`(a_s(Q^2)\left({\displaystyle \frac{a_s(Q^2)}{a_s(Q_0^2)}}\right)_2^{\frac{\lambda _+^z\lambda _{}^z}{2\beta _0}}F_1(1,1+{\displaystyle \frac{\lambda _+^z\lambda _{}^z}{2\beta _0}},2+{\displaystyle \frac{\lambda _+^z\lambda _{}^z}{2\beta _0}},{\displaystyle \frac{\beta _1}{\beta _0}}a_s(Q^2))`$
$`a_s(Q_0^2)_2F_1(1,1+{\displaystyle \frac{\lambda _+^z\lambda _{}^z}{2\beta _0}},2+{\displaystyle \frac{\lambda _+^z\lambda _{}^z}{2\beta _0}},{\displaystyle \frac{\beta _1}{\beta _0}}a_s(Q_0^2)))]`$
$`+(+).`$
Finally, the evolution operator for eq. (14) (e.g. “x-space and truncated in $`1/L`$”) reads for the non-singlet case

$`\widehat{E}_\eta ^z(Q^2,Q_0^2)`$ $`=`$ $`\left({\displaystyle \frac{L_0}{L}}\right)^{\frac{\gamma _0^z}{2\beta _0}}\mathrm{exp}[{\displaystyle \frac{1}{2\beta _0^2}}({\displaystyle \frac{1}{L_0}}{\displaystyle \frac{1}{L}})(\gamma _1^z(\eta ){\displaystyle \frac{\beta _1}{\beta _0}}\gamma _0^z)+{\displaystyle \frac{\beta _1}{2\beta _0^3}}({\displaystyle \frac{\mathrm{ln}L_0}{L_0}}{\displaystyle \frac{\mathrm{ln}L}{L}})\gamma _0^z`$ (47)
$`+{\displaystyle \frac{\beta _1}{4\beta _0^4}}\left({\displaystyle \frac{1+2\mathrm{ln}L_0}{L_0^2}}{\displaystyle \frac{1+2\mathrm{ln}L}{L^2}}\right)\gamma _1^z(\eta )`$
$`{\displaystyle \frac{\beta _1^2}{54\beta _0^6}}({\displaystyle \frac{2+6\mathrm{ln}L_0+9\mathrm{ln}^2L_0}{L_0^3}}{\displaystyle \frac{2+6\mathrm{ln}L+9\mathrm{ln}^2L}{L^3}})\gamma _1^z(\eta )],`$
where $`L_i=\mathrm{ln}Q_i^2/\mathrm{\Lambda }^2`$. In the singlet case no closed form is known and the evolution operator is obtained as a sequence of evolutions

$`\widehat{E}^z(Q^2,Q_0^2)`$ $`=`$ $`\widehat{E}^z(Q_{(n)}^2,Q_{(n1)}^2)\widehat{E}^z(Q_{(n1)}^2,Q_{(n2)}^2)\mathrm{}\widehat{E}^z(Q_{(1)}^2,Q_{(0)}^2),`$ (48)
where $`Q_{(j)}^2=Q_0^2+\left(Q^2Q_0^2\right)j/n`$. In a small interval the evolution operator for the singlet case is then approximated by

$`\widehat{E}^z(Q^2,Q_0^2)`$ $`=`$ $`\left({\displaystyle \frac{L_0}{L}}\right)^{\frac{\lambda _{}^z}{2\beta _0}}\{P_{}^z+P_{}^z[{\displaystyle \frac{1}{2\beta _0^2}}({\displaystyle \frac{1}{L_0}}{\displaystyle \frac{1}{L}})\gamma ^z+{\displaystyle \frac{\beta _1}{2\beta _0^3}}({\displaystyle \frac{\mathrm{ln}L_0}{L_0}}{\displaystyle \frac{\mathrm{ln}L}{L}})\gamma _0^z`$ (49)
$`+{\displaystyle \frac{\beta _1}{4\beta _0^4}}\left({\displaystyle \frac{1+2\mathrm{ln}L_0}{L_0^2}}{\displaystyle \frac{1+2\mathrm{ln}L}{L^2}}\right)\gamma _1^z`$
$`{\displaystyle \frac{\beta _1^2}{54\beta _0^6}}({\displaystyle \frac{2+6\mathrm{ln}L_0+9\mathrm{ln}^2L_0}{L_0^3}}{\displaystyle \frac{2+6\mathrm{ln}L+9\mathrm{ln}^2L}{L^3}})\gamma _1^z]P_{}^z`$
$`+P_{}^z\gamma _1^zP_+^z[{\displaystyle \frac{1}{\beta _0\left(2\beta _0+\lambda _+^z\lambda _{}^z\right)}}({\displaystyle \frac{1}{L_0}}{\displaystyle \frac{1}{L}}\left({\displaystyle \frac{L_0}{L}}\right)^{\delta ^z})`$
$`+{\displaystyle \frac{2\beta _1}{\beta _0^3\left(4\beta _0+\lambda _+^z\lambda _{}^z\right)}}\left({\displaystyle \frac{\mathrm{ln}L_0}{L_0^2}}{\displaystyle \frac{\mathrm{ln}L}{L^2}}\left({\displaystyle \frac{L_0}{L}}\right)^{\delta ^z}\right)`$
$`+{\displaystyle \frac{2\beta _1}{\beta _0^2\left(4\beta _0+\lambda _+^z\lambda _{}^z\right)}}\left({\displaystyle \frac{1}{L_0^2}}{\displaystyle \frac{1}{L^2}}\left({\displaystyle \frac{L_0}{L}}\right)^{\delta ^z}\right)`$
$`{\displaystyle \frac{\beta _1^2}{\beta _0^5\left(6\beta _0+\lambda _+^z\lambda _{}^z\right)}}\left({\displaystyle \frac{\mathrm{ln}^2L_0}{L_0^3}}{\displaystyle \frac{\mathrm{ln}^2L}{L^3}}\left({\displaystyle \frac{L_0}{L}}\right)^{\delta ^z}\right)`$
$`{\displaystyle \frac{2\beta _1^2}{\beta _0^4\left(6\beta _0+\lambda _+^z\lambda _{}^z\right)}}\left({\displaystyle \frac{\mathrm{ln}L_0}{L_0^3}}{\displaystyle \frac{\mathrm{ln}L}{L^3}}\left({\displaystyle \frac{L_0}{L}}\right)^{\delta ^z}\right)`$
$`{\displaystyle \frac{2\beta _1^2}{\beta _0^3\left(6\beta _0+\lambda _+^z\lambda _{}^z\right)}}({\displaystyle \frac{1}{L_0^3}}{\displaystyle \frac{1}{L^3}}\left({\displaystyle \frac{L_0}{L}}\right)^{\delta ^z})]\}+(+),`$
where $`\delta ^z=(\lambda _+^z\lambda _{}^z)/2\beta _0`$. Note that there is a typo in eq. (8.27) of ref. . The prefactor $`r_L^{\lambda _{}^z/\beta _0}`$ should read $`r_L^{+\lambda _{}^z/\beta _0}`$.