text:
# 1 Introduction

## 1 Introduction

It is known that the sunrise graph of Figure 1 in the arbitrary mass case has four master integrals;

they can be taken to be the fully scalar amplitude

$$F_0(d,m_1^2,m_2^2,m_3^2,p^2)=\frac{d^dk_1}{(2\pi )^{d2}}\frac{d^dk_2}{(2\pi )^{d2}}\frac{1}{[(pk_1)^2+m_1^2][(k_1k_2)^2+m_2^2](k_2^2+m_3^2)},$$
(1.1)
where $`d`$ is the continuous dimension, $`m_i,i=1,3`$ the masses, $`p^2`$ the square of the external momentum (all momenta being Euclidean), and the three amplitudes $`F_i`$ defined as

$$F_i(d,m_1^2,m_2^2,m_3^2,p^2)=\frac{}{m_i^2}F_0(d,m_1^2,m_2^2,m_3^2,p^2),$$
(1.2)
for $`i=1,3`$. The higher mass derivatives of the amplitudes can all be expressed, as a consequence of the integration by parts identities, in terms of the four master integrals. The four master integrals satisfy further a system of linear non homogeneous differential equations in $`p^2`$, derived in .

In the particular kinematical configuration $`p^2=m_1^2=m^2`$ and $`m_2^2=m_3^2=M^2`$, Figure 2, one has

$$F_0(d,m^2,M^2,M^2,p^2=m^2)=\frac{d^dk_1}{(2\pi )^{d2}}\frac{d^dk_2}{(2\pi )^{d2}}\frac{1}{(k_1^22pk_1)[(k_1k_2)^2+M^2](k_2^2+M^2)};$$
(1.3)
introducing the dimensionless variables $`x=m/M`$ one can define the dimensionless function $`\mathrm{\Phi }(d,x)`$ as

$$F_0(d,m^2,M^2,M^2,m^2)=M^{2d6}C^2(d)\mathrm{\Phi }(d,x),$$
(1.4)
i.e.

$$\mathrm{\Phi }(d,x)=C^2(d)M^{62d}F_0(d,M^2x^2,M^2,M^2,M^2x^2),$$
(1.5)
where $`C(d)=(4\pi )^{\frac{4d}{2}}\mathrm{\Gamma }(3d/2)`$ is an overall normalization factor, with the limiting value $`C(4)=1`$ at $`d=4`$.

It will be shown in this paper that the integration by part identities for the master integrals imply a second order differential equation in $`x`$ for $`\mathrm{\Phi }(d,x)`$; when functions and equations are systematically expanded around $`d=4`$ as a Laurent series in $`(d4)`$, one obtains a system of chained differential equations which can be easily solved recursively in closed analytical form, virtually up to any order in $`(d4)`$; one finds in this way that the coefficients of the $`(d4)`$ expansion of $`\mathrm{\Phi }(d,x)`$ are a combination of harmonic polylogarithms  of increasing weight.

The differential equations themselves with the knowledge of the behaviour at $`m=0`$ can be used to fix completely the integration constants; as a byproduct we give the value of the integral at $`m=M`$ (on mass shell equal mass case) up to the order $`(d4)^5`$ included.

## 2 The differential equations.

For deriving the differential equation for $`\mathrm{\Phi }(d,x)`$ we introduce the auxiliary functions

$`\mathrm{\Psi }_0(d,\eta ,p^2)`$ $`=`$ $`m^{62d}F_0(d,m^2,m^2\eta ,m^2\eta ,p^2),`$ (2.1)
$`\mathrm{\Psi }_1(d,\eta ,p^2)`$ $`=`$ $`m^{82d}\left[F_2(d,m^2,m^2\eta ,m^2\eta ,p^2)+F_3(d,m^2,m^2\eta ,m^2\eta ,p^2)\right].`$ (2.2)
One finds at once

$$\frac{}{\eta }\mathrm{\Psi }_0(d,\eta ,p^2)=\mathrm{\Psi }_1(d,\eta ,p^2),$$
(2.3)
from which in turn

$`{\displaystyle \frac{}{\eta }}\mathrm{\Psi }_1(d,\eta ,p^2)`$ $`=`$ $`m^{(102d)}`$ (2.4)
$`\left({\displaystyle \frac{}{m_2^2}}+{\displaystyle \frac{}{m_3^2}}\right)\left[F_2(d,m^2,m_2^2,m_3^2,p^2)+F_3(d,m^2,m_2^2,m_3^2,p^2)\right]_{m_2^2=m_3^2=m^2\eta }.`$
As already remarked, the mass derivatives of the master integrals $`F_i,i=1,3`$ can be expressed, by using the integration by parts identities, in terms of the master integrals themselves (see for instance the Appendix of ), so that the above quantity is equal to a combination of $`F_0`$, corresponding to $`\mathrm{\Psi }_0(d,\eta ,p^2)`$, of the sum $`F_2+F_3`$ (as required by the $`m_2=m_3`$ symmetry), corresponding to $`\mathrm{\Psi }_1(d,\eta ,p^2)`$, and of the fourth master integral $`F_1`$, independent of $`\mathrm{\Psi }_0`$ and $`\mathrm{\Psi }_1`$, plus an inhomogeneous part given by a few products of tadpoles $`T(d,m^2),T(d,M^2)`$, with

$$T(d,m^2)=\frac{d^dk}{(2\pi )^{d2}}\frac{1}{k^2+m^2}=C(d)\frac{m^{d2}}{(d2)(d4)}.$$
(2.5)
The coefficients of the various terms are ratios of polynomials in $`d`$, $`p^2`$ and the masses; in the limit $`p^2m^2`$ the denominators of the coefficients develop a singular $`1/(p^2+m^2)^2`$ behaviour, which is however fully compensated by corresponding zeroes in the numerators, so that the limit $`p^2=m^2`$ is finite. In that limit, further, the coefficient of the master amplitude $`F_1`$ vanishes, so that $`\frac{}{\eta }\mathrm{\Psi }_1(d,\eta ,p^2=m^2)`$, which is essentially the second $`\eta `$-derivative of $`\mathrm{\Psi }_0(d,\eta ,p^2=m^2)`$, is expressed in terms of $`\mathrm{\Psi }_1(d,\eta ,p^2=m^2)`$, its first $`\eta `$-derivative (up to a sign), and of $`\mathrm{\Psi }_0(d,\eta ,p^2=m^2)`$ itself, plus an inhomogeneous term given by the tadpoles. For taking the $`p^2=m^2`$ limit, let us introduce the function $`\mathrm{\Psi }(d,\eta )`$ through

$$F_0(d,m^2,M^2,M^2,m^2)=m^{2d6}C^2(d)\mathrm{\Psi }(d,\eta ),$$
(2.6)
i.e.

$$\mathrm{\Psi }(d,\eta )=C^2(d)m^{62d}F_0(d,m^2,m^2\eta ,m^2\eta ,m^2).$$
(2.7)
In terms of $`\mathrm{\Psi }(d,\eta )`$ the second order differential equation for $`\mathrm{\Psi }_0(d,\eta ,p^2=m^2)`$ then reads

$`\eta (1\eta ){\displaystyle \frac{^2}{\eta ^2}}\mathrm{\Psi }(d,\eta )\left[{\displaystyle \frac{1}{2}}{\displaystyle \frac{3}{2}}\eta +(12\eta )(d4)\right]{\displaystyle \frac{}{\eta }}\mathrm{\Psi }(d,\eta )`$
$`\left[1+{\displaystyle \frac{7}{4}}(d4)+{\displaystyle \frac{3}{4}}(d4)^2\right]\mathrm{\Psi }(d,\eta )`$ $`=`$ $`{\displaystyle \frac{\eta ^{(d4)/2}}{4(d4)^2}}{\displaystyle \frac{\eta ^{1+(d4)}}{8(d4)^2}}.`$ (2.8)
The equation, exact in the continuous dimension $`d`$, has been written in terms of powers of $`(d4)`$ for convenience of later use.

As $`\eta =1/x^2`$, by comparing the definitions Eq.s (1.4, 2.6) one has

$$\mathrm{\Psi }\left(d,\eta =\frac{1}{x^2}\right)=x^{62d}\mathrm{\Phi }(d,x),$$
(2.9)
and the above equation for $`\mathrm{\Psi }(d,\eta )`$ can be converted in the following equation for $`\mathrm{\Phi }(d,x)`$

$`{\displaystyle \frac{1}{4}}x^2(x^21){\displaystyle \frac{^2}{x^2}}\mathrm{\Phi }(d,x){\displaystyle \frac{1}{2}}x\left[1+(d4)x^2\right]{\displaystyle \frac{}{x}}\mathrm{\Phi }(d,x)`$
$`+\left[{\displaystyle \frac{1}{2}}(1x^2)+\left({\displaystyle \frac{3}{4}}{\displaystyle \frac{1}{2}}x^2\right)(d4)+{\displaystyle \frac{1}{4}}(d4)^2\right]\mathrm{\Phi }(d,x)`$ $`=`$ $`{\displaystyle \frac{x^{2+(d4)}}{4(d4)^2}}{\displaystyle \frac{1}{8(d4)^2}}.`$ (2.10)
## 3 The expansion in $`(d4)`$ around $`d=4`$ .

It is known  that in the $`d4`$ limit $`\mathrm{\Phi }(d,x)`$ develops a double pole in $`(d4)`$, so that its Laurent expansion in $`(d4)`$ reads

$$\mathrm{\Phi }(d,x)=\underset{n=2}{\overset{\mathrm{}}{}}(d4)^n\mathrm{\Phi }^{(n)}(x).$$
(3.1)
By inserting the above expansion in Eq. (2.10) one obtains a system of chained inhomogeneous equations in the functions $`\mathrm{\Phi }^{(n)}(x)`$. For $`n=2,1`$ they read

$`\left[{\displaystyle \frac{^2}{x^2}}+\left({\displaystyle \frac{2}{x}}+{\displaystyle \frac{1}{1x}}{\displaystyle \frac{1}{1+x}}\right){\displaystyle \frac{}{x}}{\displaystyle \frac{2}{x^2}}\right]\mathrm{\Phi }^{(2)}(x)=`$ $`{\displaystyle \frac{1}{2x^2}}{\displaystyle \frac{1}{4(1x)}}{\displaystyle \frac{1}{4(1+x)}}`$
$`\left[{\displaystyle \frac{^2}{x^2}}+\left({\displaystyle \frac{2}{x}}+{\displaystyle \frac{1}{1x}}{\displaystyle \frac{1}{1+x}}\right){\displaystyle \frac{}{x}}{\displaystyle \frac{2}{x^2}}\right]\mathrm{\Phi }^{(1)}(x)=`$ $``$ $`\left({\displaystyle \frac{1}{1x}}{\displaystyle \frac{1}{1+x}}\right){\displaystyle \frac{}{x}}\mathrm{\Phi }^{(2)}(x)`$
$`+\left({\displaystyle \frac{3}{x^2}}+{\displaystyle \frac{1}{2(1x)}}+{\displaystyle \frac{1}{2(1+x)}}\right)\mathrm{\Phi }^{(2)}(x)`$ $``$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{1x}}+{\displaystyle \frac{1}{1+x}}\right)\text{ln}x,`$ (3.2)
while for $`n0`$ one has

$`\left[{\displaystyle \frac{^2}{x^2}}+\left({\displaystyle \frac{2}{x}}+{\displaystyle \frac{1}{1x}}{\displaystyle \frac{1}{1+x}}\right){\displaystyle \frac{}{x}}{\displaystyle \frac{2}{x^2}}\right]\mathrm{\Phi }^{(n)}(x)=`$ $``$ $`\left({\displaystyle \frac{1}{1x}}{\displaystyle \frac{1}{1+x}}\right){\displaystyle \frac{}{x}}\mathrm{\Phi }^{(n1)}(x)`$
$`+\left({\displaystyle \frac{3}{x^2}}+{\displaystyle \frac{1}{2(1x)}}+{\displaystyle \frac{1}{2(1+x)}}\right)\mathrm{\Phi }^{(n1)}(x)`$ $`+`$ $`\left({\displaystyle \frac{1}{x^2}}+{\displaystyle \frac{1}{2(1x)}}+{\displaystyle \frac{1}{2(1+x)}}\right)\mathrm{\Phi }^{(n2)}(x)`$ (3.3)
$``$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{1x}}+{\displaystyle \frac{1}{1+x}}\right){\displaystyle \frac{\text{ln}^{n+2}x}{(n+2)!}}.`$
A few comments are in order. Quite in general, all the equations have the form

$$\left[\frac{^2}{x^2}+\left(\frac{2}{x}+\frac{1}{1x}\frac{1}{1+x}\right)\frac{}{x}\frac{2}{x^2}\right]\mathrm{\Phi }^{(n)}(x)=R^{(n)}(x);$$
(3.4)
the homogeneous part of the equations is always

$$\left[\frac{^2}{x^2}+\left(\frac{2}{x}+\frac{1}{1x}\frac{1}{1+x}\right)\frac{}{x}\frac{2}{x^2}\right]\varphi (x)=0,$$
(3.5)
while the inhomogeneous terms $`R^{(n)}(x)`$ vary from equation to equation and involve, besides the coefficients of the $`(d4)`$-expansion of the inhomogeneous term of Eq.(2.10) which are explicitly known since the beginning, also the coefficients $`\mathrm{\Phi }^{(n1)}(x)`$ and $`\mathrm{\Phi }^{(n2)}(x)`$ of the $`(d4)`$-expansion of the unknown function $`\mathrm{\Psi }(d,x)`$. However, in a systematic bottom-up approach to the solution starting from the smallest value of $`n`$, i.e..$`n=2`$, the inhomogeneous part can also be considered as known.

Finally, all the coefficients appearing in the equation, both in the homogeneous and the inhomogeneous part, are simple combinations of the rational factors $`1/x`$, $`1/(1x)`$ and $`1/(1+x)`$; powers of $`\text{ln}x`$ are also present, from the expansion of the original inhomogeneous term.

## 4 The behaviour of the equation for $`x0^+`$.

As $`x`$ is proportional to a mass and therefore is an essentially real and positive variable, we study Eq.(2.10) in the $`x0^+`$ limit; to that aim, we try the most general expansion

$$\mathrm{\Phi }(d,x)=\underset{i}{}x^{\alpha _i}\left(\underset{n=0}{\overset{\mathrm{}}{}}a_n^{(i)}(d)x^n\right).$$
(4.1)
By substituting the above expansion in Eq.(2.10), we find for the leading exponents $`\alpha _i`$ the four allowed values $`\alpha _1=(2d)`$, $`\alpha _2=(d3)`$, corresponding to the two independent solutions of the associated homogeneous equation and $`\alpha _3=(d2)`$, $`\alpha _4=0`$, as required by the inhomogeneous terms appearing in the r.h.s. of Eq.(2.10).

By direct inspection of Eq.(1.3) and Eq.(1.4), one sees that $`\mathrm{\Phi }(d,x)`$ is finite for $`d>2`$ in the $`x0^+`$ limit, so that the first two behaviours, with exponents $`(2d)`$ and $`(d3)`$, are ruled out (i.e. the coefficients $`a_n^{(i)}(d)`$ are all vanishing for $`i=1,2`$ and any value of $`n`$ ). The coefficients of the other terms are then completely fixed by the inhomogeneous part of the equation, and an elementary calculation gives the following explicit behaviour in the $`x0^+`$ limit

$`\mathrm{\Phi }(d,x)`$ $`=`$ $`{\displaystyle \frac{1}{2(d2)(d4)^2}}\left[{\displaystyle \frac{1}{d3}}+{\displaystyle \frac{2x^2}{d(d5)}}+x^2x^{d4}\right]`$ (4.2)
$`+`$ $`𝒪(x^3).`$
When expanding in $`(d4)`$ according to Eq.(3.1), the above equation reads

$`\mathrm{\Phi }^{(2)}(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{4}}x^2+𝒪(x^3),`$
$`\mathrm{\Phi }^{(1)}(x)`$ $`=`$ $`+{\displaystyle \frac{3}{4}}+{\displaystyle \frac{5}{16}}x^2{\displaystyle \frac{x^2}{2}}\text{ln}x+𝒪(x^3),`$
$`\mathrm{\Phi }^{(0)}(x)`$ $`=`$ $`{\displaystyle \frac{7}{8}}+{\displaystyle \frac{3}{64}}x^2+{\displaystyle \frac{x^2}{4}}(\text{ln}x\text{ln}^2x)+𝒪(x^3),`$
and so on up to any desired order in $`(d4)`$.

It is to be observed that by using the differential equation Eq.(2.10) one can transform a qualitative information, the finiteness of $`\mathrm{\Phi }(d,x)`$ at $`x=0`$, into a quantitative information, such as the actual value of the $`x0^+`$ behaviour of $`\mathrm{\Phi }(d,x)`$ given by Eq.(4.2), without carrying out any explicit loop integration in Eq.(1.3).

## 5 The solution.

The algorithm for solving Eq.(3.4) goes back to Euler. Let $`\varphi _1(x)`$, $`\varphi _2(x)`$ be two independent solutions of the associated homogeneous equation (3.5), such that their Wronskian $`W(x)`$, defined as

$$W(x)=\varphi _1(x)\varphi _2^{}(x)\varphi _1^{}(x)\varphi _2(x),$$
(5.1)
does not vanish; the solution of the second order Eq.(3.4) can then be written as

$$\mathrm{\Phi }^{(n)}(x)=\varphi _1(x)\left(\mathrm{\Phi }_1^{(n)}^x\frac{dx^{}\varphi _2(x^{})}{W(x^{})}R^{(n)}(x^{})\right)+\varphi _2(x)\left(\mathrm{\Phi }_2^{(n)}+^x\frac{dx^{}\varphi _1(x^{})}{W(x^{})}R^{(n)}(x^{})\right),$$
(5.2)
where $`\mathrm{\Phi }_1^{(n)},\mathrm{\Phi }_2^{(n)}`$ are two integration constants, to be fixed by the boundary conditions Eq.(4.2) expanded in $`(d4)`$, i.e..by Eq.s(LABEL:eq:Phiexpld).

For the actual use of Euler’s algorithm, the solutions of the homogeneous equation (3.5) are needed. It turns out that the homogeneous equation associated to Eq.(2.8) at $`d=4`$

$$\left[\eta (1\eta )\frac{^2}{\eta ^2}\left(\frac{1}{2}\frac{3}{2}\eta \right)\frac{}{\eta }1\right]\psi (\eta )=0,$$
(5.3)
is a hypergeometric equation, whose solution regular at $`\eta =0`$ is the very simple hypergeometric function

$${}_{2}{}^{}F_{1}^{}(\frac{1}{2},2,\frac{1}{2};\eta )=(1\eta )^2.$$
(5.4)
With that input and the transformation (2.9) it is easy to obtain the following solutions of the homogeneous equation (3.5):

$`\varphi _1(x)`$ $`=`$ $`{\displaystyle \frac{(1x^2)^2}{x^2}},`$
$`\varphi _2(x)`$ $`=`$ $`{\displaystyle \frac{(1x^2)^2}{x^2}}\text{ln}\left({\displaystyle \frac{1+x}{1x}}\right)2{\displaystyle \frac{1+x^2}{x}},`$ (5.5)
having the Wronskian

$$W(x)=16\frac{1x^2}{x^2},$$
(5.6)
and the $`x0`$ behaviour

$`\varphi _1(x)`$ $``$ $`{\displaystyle \frac{1}{x^2}},`$
$`\varphi _2(x)`$ $``$ $`6x.`$ (5.7)
With the above explicit homogeneous solutions one can solve recursively the system of the non homogeneous equations (3.3) for $`\mathrm{\Phi }^{(n)}(x)`$, starting from $`n=2`$, by carrying out repeatedly the integration in $`x^{}`$ of (5.2) in close analytic form; one obtains in that way a combination of Harmonic PolyLogarithms (HPLs) of increasing weight  times simple coefficients involving the rational factors $`1/x`$, $`1/(1x)`$ and $`1/(1+x)`$.

To illustrate this statement, which is the central point of this paper, let us recall for convenience of the reader the definition of the HPL’s. The HPL’s form a family of functions depending on an argument, say $`x`$, and on a set of indices, say $`a_i,i=1,..,w`$ or $`\stackrel{}{a}`$ in more compact notation, where each of the $`a_i`$ can take one of the three values $`1,0,1`$ and whose number $`w`$ is called the weight of the HPL. At weight $`w=1`$ there are 3 HPL’s, defined as

$`\text{H}(1;x)`$ $`=`$ $`\text{ln}(1x),`$
$`\text{H}(0;x)`$ $`=`$ $`\text{ln}x,`$
$`\text{H}(1;x)`$ $`=`$ $`\text{ln}(1+x),`$ (5.8)
whose derivatives can be written as

$$\frac{d}{dx}\text{H}(a;x)=\text{f}(a,x),$$
(5.9)
where the index $`a`$ can take one of the three values $`(1,0,1)`$ and the rational factors $`\text{f}(a;x)`$ are given by

$`\text{f}(1;x)`$ $`=`$ $`{\displaystyle \frac{1}{1x}},`$
$`\text{f}(0;x)`$ $`=`$ $`{\displaystyle \frac{1}{x}},`$
$`\text{f}(1;x)`$ $`=`$ $`{\displaystyle \frac{1}{1+x}}.`$ (5.10)
At weight $`w>1`$, if all the $`w`$ indices are equal to $`0`$ let us indicate them by $`\stackrel{}{0}_w`$ and define correspondingly

$$\text{H}(\stackrel{}{0}_w;x)=\frac{1}{w!}\text{ln}^wx;$$
(5.11)
in all the other cases (i.e..when the indices are not all equal to zero), let us indicate any set of $`w`$ indices by $`(a,\stackrel{}{b})`$, where $`a`$ can take one of the three values $`(1,0,1)`$ and $`\stackrel{}{b}`$ stands for the set of other $`w1`$ indices, and define correspondingly

$$\text{H}(a,\stackrel{}{b};x)=_0^x𝑑x^{}\text{f}(a;x^{})\text{H}(\stackrel{}{b};x^{}).$$
(5.12)
Note that in full generality one has

$$\frac{d}{dx}\text{H}(a,\stackrel{}{b};x)=\text{f}(a,x)\text{H}(\stackrel{}{b};x),$$
(5.13)
which can also be written as the equivalent indefinite integration formula

$$^x𝑑x^{}\text{f}(a;x^{})\text{H}(\stackrel{}{b};x^{})=A+\text{H}(a,\stackrel{}{b};x),$$
(5.14)
where $`A`$ is an integration constant.

Further, the product of two HPL’s of a same argument $`x`$ and weights $`p,q`$ can be expressed as a combination of HPL’s of argument $`x`$ and weight $`r=p+q`$, according to the product identity

$`\text{H}(\stackrel{}{p};x)\text{H}(\stackrel{}{q};x)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{r}=\stackrel{}{p}\stackrel{}{q}}{}}\text{H}(\stackrel{}{r};x),`$ (5.15)
where $`\stackrel{}{p},\stackrel{}{q}`$ stand for the $`p`$ and $`q`$ components of the indices of the two HPL’s, while $`\stackrel{}{p}\stackrel{}{q}`$ represents all possible mergers of $`\stackrel{}{p}`$ and $`\stackrel{}{q}`$ into the vector $`\stackrel{}{r}`$ with $`r`$ components, in which the relative orders of the elements of $`\stackrel{}{p}`$ and $`\stackrel{}{q}`$ are preserved. The simplest cases of the above identities are

$`\text{H}(a;x)\text{H}(b;x)`$ $`=`$ $`\text{H}(a,b;x)+\text{H}(b,a;x),`$
$`\text{H}(a;x)\text{H}(b,c;x)`$ $`=`$ $`\text{H}(a,b,c;x)+\text{H}(b,a,c;x)+\text{H}(b,c,a;x);`$ (5.16)
more complicated cases are immediately established recursively (all the above formulae can indeed easily be checked by differentiating, repeatedly when needed, with respect to $`x`$).

After this digression, let us come back to Eq.(5.2), which gives the solution as an indefinite integral. The original $`n=2`$ integrand consists of powers of $`x^{}`$, $`1/(1x^{})`$ and $`1/(1+x^{})`$, (without any loss of generality $`(1x^{})`$ and $`(1+x^{})`$ can be taken to occur only as negative powers), which are present in all the factors (see the first of Eq.s(3.2) for the definition of $`R^{(2)}(x^{})`$), and of the two HPL’s of weight 1, $`\text{H}(1;x^{})=\text{ln}(1x^{})`$ and $`\text{H}(1;x^{})=\text{ln}(1+x^{})`$, which occur in $`\varphi _2(x^{})`$, Eq.(5.5). The terms containing powers of the rational factors different from the mere inverse, i.e..different from the first powers of $`1/x^{}`$, $`1/(1x^{})`$ and $`1/(1+x^{})`$, can be integrated by parts; in so doing they give rise at most to derivatives of the occurring HPL’s, i.e..again products of the three rational factors $`x^{}`$, $`1/(1x^{})`$ and $`1/(1+x^{})`$, and, according to Eq.s(5.13, 5.9), of HPL’s of lower weight if the HPL’s in the integrand have weight greater than 1 (as it is in the general case $`n>2`$). The result can be partial fractioned, and the integration by parts repeated until all powers of the rational fractions different from the mere first power are processed. The terms with the mere inverse of the three rational fractions and HPL’s, finally, can be integrated at once according to Eq.(5.14), giving rise to HPL’s of higher weight. The integration constants are then fixed by imposing the $`x0^+`$ behaviour (LABEL:eq:Phiexpld). It is clear that the argument applies to any of the subsequent integrations as well, recursively in $`n`$.

For $`n=2`$ the explicit calculation gives $`\mathrm{\Phi }_1^{(2)}=0`$, $`\mathrm{\Phi }_2^{(2)}=0`$ and

$$\mathrm{\Phi }^{(2)}(x)=\frac{2+x^2}{8}.$$
(5.17)
Note that, in this case, all the HPL’s cancel out from the final result.

One can then proceed to $`n=1`$; the term $`\text{ln}x`$ appearing in the inhomogeneous part of the second of Eq.s(3.2) can be written as $`H(0;x)`$, see (5.8), the products of the HPL’s already present in $`R^{(1)}(x^{})`$ and the HPL’s of the $`\varphi _i(x^{})`$ can be rewritten as combinations of single HPL’s of suitable indices according to 5.15, the integration performed as explained above – and so on for the higher values of $`n`$, up to any required order. An explicit calculation gives for the integration constants $`\mathrm{\Phi }_1^{(1)}=5/64`$, $`\mathrm{\Phi }_1^{(0)}=29/256`$, $`\mathrm{\Phi }_1^{(1)}=107/1024`$, $`\mathrm{\Phi }_1^{(2)}=185/4096`$, $`\mathrm{\Phi }_1^{(3)}=1285/16384`$, $`\mathrm{\Phi }_1^{(4)}=18991/65536`$, $`\mathrm{\Phi }_1^{(5)}=164173/262144`$, while, for any $`n`$, $`\mathrm{\Phi }_2^{(n)}=0`$. The first $`\mathrm{\Phi }^{(n)}(x)`$ are

$$\mathrm{\Phi }^{(1)}(x)=\frac{3}{8}+x^2\left(\frac{5}{32}\frac{1}{4}H(0;x)\right),$$
(5.18)
$`\mathrm{\Phi }^{(0)}(x)`$ $`=`$ $`{\displaystyle \frac{3}{8}}{\displaystyle \frac{1}{8}}H(0;x)+x^2\left({\displaystyle \frac{11}{128}}+{\displaystyle \frac{5}{16}}H(0;x){\displaystyle \frac{1}{8}}H(0;x)H(0;x)\right)`$ (5.19)
$`{\displaystyle \frac{(1x^2)^2}{8x^2}}\left(H(0;x)H(1;x)H(0;x)H(1;x)H(0,1;x)+H(0,1;x)\right),`$
$`\mathrm{\Phi }^{(1)}(x)`$ $`=`$ $`{\displaystyle \frac{(1x^2)^2}{32x^2}}\times `$ (5.20)
$`\{{\displaystyle \frac{}{}}5[{\displaystyle \frac{}{}}H(0;x)H(1;x)H(0;x)H(1;x)H(0,1;x)+H(0,1;x)]`$
$`+2[{\displaystyle \frac{}{}}H(0;x)H(0;x)H(1;x)H(0;x)H(0;x)H(1;x)H(0;x)H(1;x)H(1;x)`$
$`{\displaystyle \frac{}{}}+4H(0;x)H(1;x)H(1;x)H(0;x)H(1;x)H(1;x)]`$
$`+4\left[{\displaystyle \frac{}{}}H(0,1;x)H(1;x)+2H(0,1;x)H(0;x)2H(0,1;x)H(1;x)\right]`$
$`+4\left[{\displaystyle \frac{}{}}H(0,1;x)H(1;x)2H(0,1;x)H(1;x)2H(0,1;x)H(0;x)\right]`$
$`+8[H(0,1,1;x)+H(0,1,1;x){\displaystyle \frac{3}{2}}H(0,0,1;x)+{\displaystyle \frac{3}{2}}H(0,0,1;x)`$
$`{\displaystyle \frac{1}{2}}H(0,1,1;x){\displaystyle \frac{1}{2}}H(0,1,1;x)]{\displaystyle \frac{}{}}\}`$
$`+`$ $`{\displaystyle \frac{(1+x)^2}{x}}\left[{\displaystyle \frac{}{}}H(0,1;x)H(0;x)H(1;x)\right]`$
$`+`$ $`{\displaystyle \frac{(1x)^2}{x}}\left[{\displaystyle \frac{}{}}H(0,1;x)H(0;x)H(1;x)\right]`$
$``$ $`{\displaystyle \frac{x^2}{8}}\left[{\displaystyle \frac{55}{64}}+{\displaystyle \frac{11}{8}}H(0;x){\displaystyle \frac{5}{4}}H(0;x)H(0;x)+{\displaystyle \frac{1}{3}}H(0;x)H(0;x)H(0;x)\right]`$
$`+`$ $`{\displaystyle \frac{15}{64}}+{\displaystyle \frac{13}{32}}H(0;x){\displaystyle \frac{1}{16}}H(0;x)H(0;x).`$
Eq.s(5.17,5.18) are of course in agreement with the general mass case (see for instance ,); Eq.(5.19) reproduces Eq.(25) of , while Eq.(5.20) is new in the literature. As matter of fact we evaluated explicitly, by means of an integration routine written in FORM , all the $`\mathrm{\Phi }^{(n)}(x)`$ up to $`n=5`$ included, which is found to involve HPL’s of weight $`\mathrm{𝟕}`$. The explicit expressions obtained are however too long to be listed here. It is found that, in general, $`\mathrm{\Phi }^{(n)}(x)`$ involves HPL’s of weight $`(n+2)`$.

The values of the HPL’s of argument equal to 1 are known ; we can then take our explicit analytic results for the $`\mathrm{\Phi }^{(n)}(x)`$ up to $`n=5`$ and evaluate them at $`x=1`$ by using , obtaining the analytical on-shell value of the sunrise amplitude Eq. 1.3 in the equal mass case limit up to the 5th order in the $`(d4)`$ expansion. At $`x=1`$ the terms of highest weight disappear from $`\mathrm{\Phi }^{(n)}(x)`$, due to the general structure of the solution, Eq.(5.5,5.2), and the coefficient of $`(d4)^n`$ involves constants up to weight $`n+1`$ only. The result reads

$`\mathrm{\Phi }(d,x=1)`$ $`=`$ $`{\displaystyle \frac{3}{8}}{\displaystyle \frac{1}{(d4)^2}}+{\displaystyle \frac{17}{32}}{\displaystyle \frac{1}{(d4)}}{\displaystyle \frac{59}{128}}+(d4)\left({\displaystyle \frac{65}{512}}+{\displaystyle \frac{1}{24}}\pi ^2\right)`$ (5.21)
$`+(d4)^2\left({\displaystyle \frac{1117}{2048}}{\displaystyle \frac{13}{96}}\pi ^2+{\displaystyle \frac{1}{8}}\pi ^2\text{ln}2{\displaystyle \frac{7}{16}}\zeta (3)\right)`$
$`+(d4)^3\left({\displaystyle \frac{13783}{8192}}+{\displaystyle \frac{115}{384}}\pi ^2{\displaystyle \frac{13}{32}}\pi ^2\text{ln}2+{\displaystyle \frac{91}{64}}\zeta (3)+{\displaystyle \frac{1}{8}}\pi ^2\text{ln}^22{\displaystyle \frac{31}{2880}}\pi ^4+{\displaystyle \frac{1}{16}}\text{ln}^42+{\displaystyle \frac{3}{2}}a_4\right)`$
$`+(d4)^4({\displaystyle \frac{114181}{32768}}{\displaystyle \frac{865}{1536}}\pi ^2+{\displaystyle \frac{115}{128}}\pi ^2\text{ln}2{\displaystyle \frac{805}{256}}\zeta (3){\displaystyle \frac{13}{32}}\pi ^2\text{ln}^22+{\displaystyle \frac{403}{11520}}\pi ^4{\displaystyle \frac{13}{64}}\text{ln}^42`$
$`{\displaystyle \frac{39}{8}}a_4{\displaystyle \frac{31}{960}}\pi ^4\text{ln}2+{\displaystyle \frac{1}{8}}\pi ^2\text{ln}^32+{\displaystyle \frac{5}{96}}\pi ^2\zeta (3)+{\displaystyle \frac{3}{80}}\text{ln}^52{\displaystyle \frac{9}{2}}a_5+{\displaystyle \frac{465}{128}}\zeta (5))`$
$`+(d4)^5({\displaystyle \frac{820495}{131072}}+{\displaystyle \frac{5971}{6144}}\pi ^2{\displaystyle \frac{865}{512}}\pi ^2\text{ln}2+{\displaystyle \frac{6055}{1024}}\zeta (3)+{\displaystyle \frac{115}{128}}\pi ^2\text{ln}^22{\displaystyle \frac{713}{9216}}\pi ^4+{\displaystyle \frac{115}{256}}\text{ln}^42`$
$`+{\displaystyle \frac{345}{32}}a_4{\displaystyle \frac{65}{384}}\pi ^2\zeta (3){\displaystyle \frac{13}{32}}\pi ^2\text{ln}^32+{\displaystyle \frac{403}{3840}}\pi ^4\text{ln}2{\displaystyle \frac{39}{320}}\text{ln}^52{\displaystyle \frac{6045}{512}}\zeta (5)+{\displaystyle \frac{117}{8}}a_5`$
$`{\displaystyle \frac{9}{160}}\pi ^4\text{ln}^22+{\displaystyle \frac{1}{64}}\pi ^2\text{ln}^42+{\displaystyle \frac{79}{34560}}\pi ^6+{\displaystyle \frac{15}{8}}\zeta (3)\text{ln}^32{\displaystyle \frac{25}{32}}\zeta (3)\pi ^2\text{ln}2+{\displaystyle \frac{595}{128}}\zeta ^2(3)`$
$`+{\displaystyle \frac{45}{4}}\zeta (5)\text{ln}2{\displaystyle \frac{45}{4}}a_5\text{ln}2+{\displaystyle \frac{21}{320}}\text{ln}^629a_6+{\displaystyle \frac{45}{4}}b_6)`$
$`+𝒪((d4)^6),`$
where we have expressed the constant $`\mathrm{𝚜𝟼}`$ appearing in  as

$`\mathrm{𝚜𝟼}`$ $`=`$ $`{\displaystyle \frac{1}{6}}\pi ^2\zeta (3)\text{ln}2{\displaystyle \frac{1}{72}}\pi ^2\text{ln}^42{\displaystyle \frac{1}{720}}\pi ^4\text{ln}^22+{\displaystyle \frac{1}{540}}\pi ^6+2\zeta (5)\text{ln}2`$ (5.22)
$`2a_5\text{ln}2+{\displaystyle \frac{1}{3}}\zeta (3)\text{ln}^32+{\displaystyle \frac{1}{120}}\text{ln}^62+{\displaystyle \frac{5}{4}}\zeta ^2(3)4a_6+2b_6,`$
and the definitions of all the previous constants in terms of HPL’s or Nielsen’s polylogarithms (NPl) is given in Table 1, together with their numerical values.

Our result, Eq.(5.21), agrees with previous results in the literature (apart from obvious normalization factors), such as , which goes up the $`(d4)^3`$ of present work, and  which arrives at $`(d4)^4`$.

Numerically, Eq.(5.21) gives

$`\mathrm{\Phi }(d,x=1)`$ $`=`$ $`0.37500000000000000000{\displaystyle \frac{1}{(d4)^2}}+0.53125000000000000000{\displaystyle \frac{1}{(d4)}}`$ (5.23)
$`0.46093750000000000000+0.53818664171204166667(d4)`$
$`0.46186261021407291667(d4)^2+0.53810855601624843750(d4)^3`$
$`0.46190033063946289062(d4)^4+0.53809670843016515942(d4)^5`$
$`+𝒪((d4)^6),`$
in perfect agreement with the numerical results of , which moreover gives two further terms in the expansion.

As a last curiosity, let us observe that for $`d4`$ Eq.(5.23) can be rewritten as

$`\mathrm{\Phi }(d,x=1)`$ $`=`$ $`{\displaystyle \frac{1}{(d5)}}(+0.37500000000000000000{\displaystyle \frac{1}{(d4)^2}}0.90625000000000000000{\displaystyle \frac{1}{(d4)}}`$ (5.24)
$`+0.992187500000000000000.99912414171204166667(d4)`$
$`+1.0000492519261145833(d4)^20.99997116623032135417(d4)^3`$
$`+1.0000088866557113281(d4)^40.99999703906962805005(d4)^5`$
$`{\displaystyle \frac{}{}}+𝒪((d4)^6)),`$
showing that the coefficients of the powers of $`(d4)^n`$ approach quickly $`(1)^n`$ for increasing $`n`$.

## 6 Conclusions.

We have shown that the two-loop self-mass sunrise amplitude, in the usual $`d`$-continuous dimensional regularization and particular kinematical configuration corresponding to two equal masses $`M`$ and the external invariant equal to the square of the third mass $`m`$, satisfies a second order differential equation in the dimensionless ratio $`x=m/M`$. The equation can be expanded in Laurent series at $`d=4`$, and solved in closed analytical form. As a results, all the $`n`$-th coefficient of the Laurent expansion can be expressed in terms of harmonic polylogarithms of argument $`x`$ and weight up to $`n+2`$. As a by product, we give the explicit analytic expressions of the value of the amplitude at $`x=1`$, corresponding to the on-mass-shell sunrise amplitude in the equal mass case, up to the $`(d4)^5`$ term included.

Acknowledgments

PM was supported in part by a grant of INFN (bando 7798/99) and in part by a grant of the ”Marie Curie Training Site” at Karlsruhe University.


format: text

text:
# Charmed hadron physics in quenched anisotropic lattice QCDPoster presented by H. Matsufuru.

## 1 Introduction

Recent experimental progress in flavor physics to look for the effect of new physics strongly requires precise theoretical predictions from the standard model. However, model independent calculation of hadronic matrix elements is difficult because of nonperturbative nature of QCD. The lattice QCD simulation is one of most promising approaches in which the systematic uncertainties can be reduced systematically . The ultimate goal of this work is to construct a framework for lattice calculations of hadronic matrix elements in a few percent systematic accuracy, as required by the experiments in progress.

In a lattice calculation of heavy quarks such as the charm and bottom, one need to either avoid or control the large discretization error of $`O((m_Qa)^n)`$. Although the effective theoretical approaches  successfully describe the matrix elements with heavy quarks within 10% accuracy, further improvement of accuracy is difficult due to perturbative error. The nonperturbative renormalization method  are not directly applicable to a determination of the coefficient of $`O(a)`$ improving term in the region of $`m_Qa^1`$ because of a strong mass dependence. Therefore a new framework is desired for which a systematic nonperturbative improvement can be performed.

The anisotropic lattice, on which the temporal lattice spacing $`a_\tau `$ is finer than the spatial one $`a_\sigma `$, is a candidate of such a framework. A high temporal lattice cutoff may enable a relativistic treatment of heavy quark while keeping computational requirement tractable. The $`O(a)`$ improvement coefficients determined in the light quark region in good precision may directly be available for a heavy quark with $`m_Qa_\tau ^1`$. Whether these observations practically hold should be examined numerically, as well as in the perturbation theory.

In this study, we examine the applicability of anisotropic lattice approach numerically by computing the charmed meson masses and decay constants. Although presently achieved accuracy is far from our goal, the result is encouraging for further development in this direction. This paper shows our preliminary result without technical details and evaluation of the systematic uncertainties, which are fully discussed in other publications .

## 2 Anisotropic lattice quark action

The quark action we adopt has the following structure in the hopping parameter form :

$`S_F`$ $`=`$ $`{\displaystyle \underset{x,y}{}}\overline{\psi }\left(x\right)\{\delta _{x,y}\kappa _\tau [(1\gamma _4)U_4\left(x\right)\delta _{x+\widehat{4},y}`$ (1)
$`+(1+\gamma _4)U_4^{}(x\widehat{4})\delta _{x\widehat{4},y}]`$
$`\kappa _\sigma _i\left[\left(r\gamma _i\right)U_i\left(x\right)\delta _{x+\widehat{i},y}\left(r+\gamma _i\right)U_i^{}\left(x\widehat{i}\right)\delta _{x\widehat{i},y}\right]`$
$`\kappa _\sigma c_E_i\sigma _{4i}F_{4i}\left(x\right)\delta _{x,y}`$
$`r\kappa _\sigma c_B_{i>j}\sigma _{ij}F_{ij}\left(x\right)\delta _{x,y},\left\}\psi \right(y),`$
where $`\kappa _\sigma `$ and $`\kappa _\tau `$ are the spatial and temporal hopping parameters, and are related to the bare quark mass and bare anisotropy as given below. $`r`$ is the spatial Wilson parameter, and $`c_E`$ and $`c_B`$ are the coefficients of clover terms which remove $`O(a)`$ errors. Although the explicit Lorentz symmetry is removed due to the anisotropy, in principle it can be restored for physical observables up to errors of $`O(a^2)`$ by proper tuning of $`\kappa _\sigma /\kappa _\tau `$, $`r`$, $`c_E`$ and $`c_B`$ for a given $`\kappa _\sigma `$. The action is constructed in accord with the Fermilab approach  and hence applicable to an arbitrary quark mass. However, a mass dependent tuning of parameters is difficult beyond perturbation theory. This may be circumvented by taking $`a_\tau ^1m_Q`$, with which the mass dependence of parameters are expected to be small, and it suffices with those values in the chiral limit.

In present study, we vary only two parameters $`\kappa _\sigma `$ and $`\kappa _\tau `$ with fixed other parameters. We set the Wilson parameter as $`r=1/\xi `$ and the clover coefficients as the tadpole-improved tree-level values, $`c_E=1/u_\sigma u_\tau ^2`$, and $`c_B=1/u_\sigma ^3`$. The tadpole improvement  is achieved by rescaling the link variable as $`U_i(x)U_i(x)/u_\sigma `$ and $`U_4(x)U_4(x)/u_\tau `$, with the mean-field values of the spatial and temporal link variables, $`u_\sigma `$ and $`u_\tau `$, respectively. Instead of $`\kappa _\sigma `$ and $`\kappa _\tau `$, we introduce $`\kappa `$ and $`\gamma _F`$ as

$`{\displaystyle \frac{1}{\kappa }}`$ $``$ $`{\displaystyle \frac{1}{\kappa _\sigma u_\sigma }}2(\gamma _F+3r4)=2(m_0\gamma _F+4),`$
$`\gamma _F`$ $``$ $`\kappa _\tau u_\tau /\kappa _\sigma u_\sigma .`$ (2)
The former plays the same role as on the isotropic lattice, and the latter corresponds to the bare anisotropy.

On an anisotropic lattice, one must tune the parameters so that the anisotropy of quark field, $`\xi _F`$, equals that of the gauge field, $`\xi _G`$:

$$\xi _F(\beta ,\gamma _G;\kappa ,\gamma _F)=\xi _G(\beta ,\gamma _G;\kappa ,\gamma _F)=\xi .$$
(3)
Although $`\xi _G`$ and $`\xi _F`$ are in general functions of both of gauge parameters ($`\beta `$, $`\gamma _G`$) and quark parameters ($`\kappa `$, $`\gamma _F`$), on a quenched lattice one can determine $`\xi _G=\xi `$ independently of $`\kappa `$ and $`\gamma _F`$, and then tune $`\gamma _F`$ so that a certain observable satisfies the condition (3). In this work, we define $`\xi _F`$ through the relativistic dispersion relation of meson,

$$E^2(𝒑)=m^2+𝒑^2/\xi _F^2+O[(𝒑^2)^2],$$
(4)
for calibration. In the above expression, the energy and mass $`E`$ and $`m`$ are in temporal lattice units while the momentum $`𝒑`$ is in spatial lattice units. $`\xi _F`$ converts the momentum into that in temporal lattice units.

## 3 Numerical simulation

Numerical simulations are performed on two quenched lattices of sizes $`16^3\times 128`$ and $`20^3\times 160`$ with the Wilson plaquette action at $`\beta =5.95`$ and $`6.10`$, respectively, with the renormalized anisotropy $`\xi =4`$ . The values of bare anisotropy are chosen according to a numerical result performed in one percent accuracy in . The lattice scales are set by $`K^{}`$ meson mass, and result in $`a_\sigma ^1=1.525(27)`$ GeV and $`1.817(22)`$ GeV for $`\beta =5.95`$ and $`6.10`$, respectively. These values deviate to about 10% from scales determined with other physical quantity, say hadronic radius $`r_0`$, and represent the effect in neglecting dynamical quarks. The mean field values $`u_\tau `$ and $`u_\sigma `$ are obtained in the Landau gauge .

The calibration was done in the quark mass region below around charm quark mass in . It was found that the quark mass dependence of bare anisotropy $`\gamma _F^{}`$, with which $`\xi _F=\xi `$ holds, is actually small, and $`\gamma _F^{}`$ is well fitted to a linear form in $`m_q^2`$,

$$\frac{1}{\gamma _F^{}}(m_q)=\zeta _0+\zeta _2m_q^2,m_q=\frac{1}{2\xi }\left(\frac{1}{\kappa }\frac{1}{\kappa _c}\right).$$
(5)
The fit results in the values ($`\zeta _0,\zeta _2,\kappa _c`$) = (0.2490(8), 0.189(15), 0.12592(8)) at $`\beta =5.95`$, and (0.2479(9), 0.143(14), 0.12558(9)) at $`\beta =6.10`$. These values are determined in 1% statistical accuracy, while in the chiral limit additional 1% error exists due to the form of extrapolation. It was also shown that the systematic discretization errors decrease toward the continuum limit. The light hadron spectra with the obtained parameters are consistent with previous works on isotropic lattices.

For heavy quark, we use four values of $`\kappa `$ covering the physical charm quark mass with $`\gamma _F`$ according to Eq. (5). As the light quark, we use three values of $`\kappa `$ with $`\gamma _F`$ at the massless limit, which correspond to the lightest three values used in  for light hadron spectroscopy. The meson correlator is calculated with the local meson operator. For the decay constant, we focus on the pseudoscalar mesons. While one needs the matching constants to convert the lattice result into the continuum theory, we only perform the tadpole-improved tree-level matching in this paper.

We first extrapolate the result to the chiral limit, linearly in light pseudoscalar meson mass squared. The physical ($`u`$, $`d`$) and $`s`$ quark masses are defined with the massless limit and the $`K`$ meson mass. Then they are fitted to the quadratic form in $`1/m_H`$, $`m_H`$ the PS heavy-light meson mass, and interpolated to the physical $`m_D`$ meson mass. The decay constant of heavy-light meson is extrapolated and interpolated in the form of $`f_H\sqrt{m_H}`$. The result after chiral extrapolation is shown in Figure 1. The gross feature of $`m_H`$ dependence is consistent with previous works . The masses and decay constants for physical quark masses are listed in Table 1.

The hyperfine splittings are consistent with previous works, while 35% less than the experimental values. On the other hand, the splitting $`m_{D_s}m_D`$ is close to the experimental value. The $`\beta `$ dependence of decay constants are rather large, and considered as of the discretization error and renormalization effect, as well as the uncertainty in setting the scale on quenched lattice. Taking ratios of decay constants, $`\beta `$ dependence is largely cancelled. The result is consistent with the previous lattice works . For more quantitative discussion, we need to discuss more carefully on the systematic uncertainties due to the anisotropy, and to include the matching of current with continuum theory. As our conclusion, the result of numerical simulation is encouraging to pursuit more quantitative analysis in this direction.


text:
# Total-energy-based structure prediction for d(AlNiCo)

## I Introduction

This paper reports on early results of a program to determine the structure of a quasicrystal – specifically $`d`$(AlNiCo) – strictly from energetics, without use of diffraction data. Almost all of the work was performed by M.M. and M.W., with their coauthors, and is reported in Ref. , but the present account is tilted towards the particular interests of C.L.H.

We have focused on the thermodynamically stable and highly ordered “Basic Nickel” subphase of decagonal $`d`$(AlNiCo) , of composition approximately Al0.70Ni0.21Co0.09. Our ultimate objective is to predict the structure on the basis of total energy; though the gross aspects of the structure are clear from diffraction, there are many ambiguous or controversial details, and there is no understanding for the quasicrystal’s special stability and simplicity at this composition. In particular, we wish to identify the sites of the transition metal (“TM”) components Ni and Co, which are indistinguishible to (ordinary) X-rays or electrons .

In a sense, the most difficult problem of crystal chemistry is to predict which structure a given composition will adopt, even if one has an exact and tractable Hamiltonian for the total energy (e.g. by pair potentials) – for this necessarily involves a comparison with an infinity of possible structures. The only irreproachable method is a mathematical proof, which is feasible in rare cases with short-range interactions, e.g. close-packing of hard spheres, or the two-dimensional “binary-tiling” quasicrystal (a toy model) . A naive, brute-force approach would be to cool (in simulation) from the melt, and see which structures emerge; but this is prone to fail for a complicated material, since the accessible times are so short that it will get stuck in a glassy disordered configuration. Even when a quasicrystal emerges, as with Dzugutov’s toy potentials , the system may have found a merely metastable ordered state for kinetic reasons: indeed, the stable phase in this case turned out to be a simple bcc packing .

## II Decoration models

We represent the quasicrystal structure as a decoration of disjoint tiles. There are several reasons for doing so, rather than as a irrational cut through a five-dimensional hypercrystal. Based in real space, this representation is easier to visualize and somewhat more tractable, technically, than the hypercrystal approach. For example, shifts of the atoms from ideal (tiling-vertex) sites are parametrized by just a finite set of real numbers, like the coordinates in a crystal unit cell. A decoration description – particularly in the decagonal case – naturally lends itself to a hierarchy of supertilings, which are presumably involved in the modulations which distinguish some of the subphases in the $`d`$(AlNiCo) phase diagram . Furthermore, most hyperspace structures which have been defined using discrete acceptance domains, e.g. Ref. , can be translated into tiling-decoration language.

Decoration models are uniquely suited to a correct modeling of chemical disorder (including vacancies) in a quasicrystal. Of course, a hyperspace fit to Bragg data can incorporate mixed occupation of a site type, but it is highly implausible (from the viewpoint of the structural energy) that the occupations of neighboring sites vary independently of each other. If the occupations are strongly correlated, then the mean structure may differ significantly from any particular real structure; this matters e.g. because the electronic structure is quite sensitive to the TM-TM contacts in $`d`$(AlNiCo) . The tiling decoration is a very convenient way to account for such occupation correlations, by ascribing the randomness to the tiles and not the decoration. The ultimate criterion of a decorated model – fixing e.g. the appropriate size of tile – is that the ensemble of allowed tilings corresponds one-to-one with the ensemble of low-energy arrangements in an atomistic model.

A decoration description replaces the independent real coordinates of many atoms by tiling degrees of freedom, which are discrete and many times fewer. When the decoration is deterministic and disjoint (each atom bound to a particular tile or tile-related geometric object), then the atoms’ interaction energies can be expressed as a function solely of the tile configuration, called the “tile Hamiltonian” $`_{\mathrm{tile}}`$. This facilitates faster Monte Carlo simulations, on the comparatively few tile degrees of freedom, in order (i) to reveal supertilings in the process of discovering the atomic structure, as in this work; (ii) to predict diffuse scattering, or (iii) eventually to measure the phason elastic constants (which govern the diffuse wings surrounding Bragg peaks.) Armed with a reliable $`_{\mathrm{tile}}`$, one can also simulate three-dimensional samples to obtain quantitative estimates of the perfection of the long-range order and the strength of diffuse diffraction.

This approach may also be the easiest way to resolve, for a particular stable quasicrystal, whether or not it is entropically stabilized. It is stabilized by energy favoring a quasiperiodic state if $`_{\mathrm{tile}}`$ happens to implement the Penrose ‘matching rules’. That almost occurs in $`d`$(AlCuCo), as modeled with potentials like those used in the present work .

Our challenge, then, is to find the decoration rule and the tile-packing constraints without any bias, apart from the above-mentioned assumptions about the lattice constants. A caveat to keep in mind is that the results will probably be quite sensitive to the exact composition and number density.

## III Inputs: potentials and constraints

The main input to our calculations is a set of atom-atom pair potentials $`V_{ij}(r)`$ for species $`i`$ and $`j`$, six distinct functions in the case of a ternary alloy. They are given by Moriarty’s “generalized pseudopotential theory” (GPT) , a systematic expansion of the total energy as a sum $`E_0+E_1+E_2+\mathrm{}`$, where the $`E_n`$ term depends on $`n`$ atoms, but we use only the $`n=2`$ terms. The first-neighbor well of the raw $`V_{\mathrm{TMTM}}`$ is unphysically deep (which would be canceled by $`E_3`$ in the systematic theory). Therefore our $`V_{\mathrm{TMTM}}`$ potentials were empirically modified, by fitting a short-range repulsive correction so as to match all the forces in a full density-functional calculation on a small (50 atom) approximant of $`d`$(AlNiCo) . There is strong support that our $`V_{ij}(r)`$ have the correct $`r`$ dependence and relative strength, since they predict the Al-Co and Al-Ni phase diagrams pretty well  (as a function of concentration at $`T=0`$); however the ternary Al-Ni-Co phase diagram was not attempted. On the basis of the simulated melting temperature and phonon spectra, we do suspect that our potentials are $``$30% stronger than the reality.

Like the similar Al-TM potentials of Phillips, Zou, and Carlson, used earlier in Refs.  and , the GPT potentials depend implicitly on the net valence electron density, and exhibit strong “Friedel oscillations” as a function of distance (tails decaying as $`\mathrm{cos}(2k_Fr+\delta )/r^3`$). For example, $`V_{\mathrm{AlTM}}(r)`$ and $`V_{\mathrm{TMTM}}(r)`$ show four prominent minima in the range $`r8.5`$ Å; such minima are quite important in deciding the structure .

Notice that V$`{}_{\mathrm{AlAl}}{}^{}(r)`$ is practically zero after its hard core radius $`2.9`$ Å : there is not even a nearest-neighbor well, just a shoulder on the tail of the hard-core. The VAlTM potentials have very deep nearest-neighbor wells: this creates the illusion of a TM-TM repulsion, but really TM-TM neighbors are avoided mainly to make more room for Al in the TM coordination shells. However, at roughly 30% TM content, some TM-TM contacts are unavoidable; they are all Ni-Ni, since $`V_{\mathrm{AlCo}}`$ has a deeper well than VAlNi. Finally, the $`V_{\mathrm{TMTM}}`$ potentials have their deepest well at second-neighbor ($`4.2`$Å) distances. The consequence of all this is that the TM atoms form a rather rigid and somewhat uniformly spaced network , while Al atoms move rather freely to follow the potential wells or troughs created by the TM arrangement. (See Sec. V for more discussion).

Our use of these potentials imposes some limitations. We omitted the $`n=0`$ term, which depends only on the valence electron density, and is the largest contribution to a metal’s total energy. Hence we can make valid comparisons only between structures with (practically) the same valence electron density. We cannot meaningfully predict the lattice constants. Nor can we even expect the quasicrystal phase to be globally stable according to our potentials, since they might spuriously favor phase-separation.

In view of the above-mentioned pitfalls for a brute-force approach, we adopted as a second input the experimental lattice constants : the stacking period is $`c=4.08`$ Å and the quasilattice constant in the decagonal plane is $`a_0=2.45`$ Å. We only seek the lowest energy among arrangements on this framework. This highly constraining assumption still permits a vast ensemble of possible structures.

## IV Methods and results

Our methods are a mix of tile-flip Monte Carlo, atom-hopping Monte Carlo, relaxation of atomic positions, and molecular dynamics. Our procedure is first to discover the favorable low-energy motifs through Monte Carlo annealing, then to remove unnecessary degrees of freedom, and repeat, producing successively more constrained models. Having fewer degrees of freedom the latter are much faster to simulate at low temperatures. At the end, we can investigate the effects of letting atoms depart from ideal sites.

The initial stage of our exploration stacks two independent small (edge $`a_0`$) Penrose rhombic tilings, in a vertical space of one lattice constant $`c=4.08`$ Å and the only allowed atom sites are on vertices. A manageable size was 50 atoms on 72 candidate sites, initially on a good Penrose-tiling approximant; periodic boundary conditions are used in all directions (and at all stages of exploration). The atoms – initially chosen to approximate the experimental density and composition of the Ni-rich “Basic Ni” phase – hop as a lattice gas on these sites. This allows sufficient freedom for the atoms, if they “want”, to adopt any of the decagonal structure models – with stacking period $`c=4`$ Å – that were hypothesized at one time or another . The Monte Carlo moves permit swaps of the species between two nearby sites (“vacant” is a special case of species!), as well as “tile flips” which reshuffle the three rhombi in a fat or thin hexagon in the same layer.

After this model is slowly cooled to zero temperature, if the initial composition was rightly chosen, one obtains a one-layer Hexagon-Boat-Star (HBS) tiling of edge length $`a_0`$, with a two-layer decoration in which the allowed sites lie over vertices of the Penrose rhombi (into which HBS tiles may be decomposed) or an additional site in each Fat rhombus. This decoration places Al atoms over the HBS vertices (even and odd vertices alternate between even and odd layers). Isolated Co atoms sit over each Boat and Star tile center, ringed by possible Al sites; in the Star at most two Al can be present, out of five ideal sites in this ring. Particularly characteristic are the vertical zigzag chains, with Ni in each layer (appearing as NiNi pairs in Fig. 1), over the interior of every Hexagon tile. However, NiNi may be replaced by CoAl in some places. Decorating a tiling with the same H/B/S ratios as the quasiperiodic Penrose tiling yields an ideal composition Al0.700Ni0.207Co0.093 and atomic volume 14.16 Å3. This composition coincides with the experimental Basic Ni phase, but the number density of atoms is at least 5% greater than in experiment.

For the next stage of modeling, to discover larger-scale regularities, the small HBS tiles which emerged from the initial stage are elevated to fundamental objects, either with a rigid decoration or with some atoms fixed and others forming a lattice gas. The allowed tile flips are reshufflings of fat hexagons of the underlying rhombus tiling, provided the result is a valid HBS tiling and conserves the atom content. Additionally the Al pair inside the Star tile can rotate among five allowed orientations. From this simulation, it emerges that the Hexagon tiles containing Ni chains only touch tip-to-tip, so that the angles relating them are multiples of $`72^{}`$. In fact, (see Fig. 1) the long axes of the Hexagons form edges of an HBS supertiling with an inflated edge length $`\tau ^2a_0`$, where $`\tau (1+\sqrt{5})/2`$.

## V The 20 Å decagon cluster

We now compare our results with well-known structure models, organizing the discussion around the famous 20 Å diameter decagon cluster, which is prominent in electron-microscope images and in most diffraction refinements as well. (The ideal decagon diameter is actually $`2\tau ^3a_0`$.) Such decagons indeed appear in our structures ; Fig. 1 is centered on one of them. Z-contrast images, in which the intensity is a direct projection of atoms weighted by squared atomic number, reveal the TM positions ; the TM (and most Al) positions proposed by Ref.  are practically the same as in our version of the cluster.

For a more detailed comparison, we focus on three recent structural studies . All of these, and our simulation, agree on the following details of the projection: (1) the outer decagon (edge length $`\tau ^2a_0`$) has Al on each vertex and a pair of TM atoms on each edge, which we identify as Ni. (2) A middle decagon (edge $`\tau a_0`$) has TM atoms on each vertex, which we identify as Co. (3) An inner decagon (edge $`a_0`$) has Al on each vertex. (4) In the center, the 10-fold symmetry is broken and a sort of isosceles triangle is observed, with one TM (we say Co) at the unique corner and pairs of TM (we say NiNi) at the base corners. The cluster is evident in Ref.  on the right side of their Fig. 7, as a combination of a Boat tile + 2 Hexagon tiles, just as in our Fig. 1.

Detail (4) is somewhat controversial, since some decagon images have non-triangular centers. Indeed, the density maps from the refinement of Ref.  show six strongly TM sites at the center of many (but not all) of the 20 Å decagons. We ascribe this to stacking flips between one layer and the next layer (see Sec. VI), seen in projection, as in the right-side decagon of Fig. 2.

The simulations we described up till now implicitly assumed a strict $`c=4.08`$ Å periodicity in the stacking direction, neglecting (like so much other modeling) the fact that decagonal quasicrystals are three-dimensional. To obtain the configuration in Fig. 2, we increased the periodicity to $`2c`$, i.e. two HBS layers, which initially were the same tiling. First a tile-flip was made in one of the two HBS layers, then the tilings were decorated, annealed by MD, and relaxed to an energy minimum.

In a fifth, controversial detail of the decagons, our model initially disagreed with experiments in which every edge of the middle decagon shows a pair of Al atoms. These sites form vertical zigzag chains, but they cannot all be occupied since they are separated by essentially the interlayer spacing $`c/2=2.08`$ Å. In our fixed-decoration model, then, these sites are occupied in only one atom layer, with separation $`c`$, so only one atom is visible on each edge. We found that, after molecular dynamics (MD) at 1000K followed by relaxation, some Al neighbors of the Co atoms shifted to the middle-decagon edges, which now present Al doublets in projection, in agreement with experiment.

To further understand what happens with these atoms, the time-average of the Al positions during the MD simulation is shown in Fig. 3. Some 40% of Al atoms are rather delocalized, and would need a highly anisotropic Debye-Waller factor in a crystallographic fit.  The vertical projection (Fig. 3(a) shows, consistent with our remarks on the potentials (Sec. III), that each Co atom in the middle decagon is surrounded by a potential trough in which Al atoms appear almost free to roll like ball-bearings .

A slice along a vertical plane further clarifies the Al behavior: a second type of trough extends vertically, with a zigzag shape, and in fact connects with the circular troughs. In our simulations with a cell $`2c`$ in the vertical direction, we found three Al atoms appearing in each zigzag trough (per $`2c`$, i.e. per 4 atom layers). Notice that the $`z`$ displacement of two of these atoms makes the layer puckered, as is already known from diffraction , and implies a local cell-doubling along the $`c`$ axis, similar to the $`\mathrm{Al}_{13}\mathrm{Fe}_4`$ or $`\mathrm{Al}_{13}\mathrm{Co}_4`$ decagonal approximants.

The symmetry-breaking of the cluster interior was recently predicted from energies by a full ab-initio calculation . They find, from a quite different starting point than ours, the identical arrangement of the 5 Al + 5 TM atoms (per atomic bilayer) found at the cluster’s center. However, they rejected the possibility of Al doublets on the middle decagon edges without trying $`z`$ relaxations, nor did they address Co/Ni ordering. We find it striking that our mere pair potentials sufficed to obtain the experimentally indicated form of symmetry breaking in the cluster center, while using no experimental input whatsoever regarding the positioning of Al relative to TM atoms.

## VI Tile Hamiltonian and three-dimensional stacking

To model a large, three-dimensional sample using a tile Hamiltonian, we cannot demand perfect periodicity in the $`z`$ direction. (One reason is that the “entropic stabilization” explanation of quasicrystal order requires an extensive entropy.) So, we represent each bilayer of atoms by a distinct HBS tiling, and constrain adjacent tilings in the stack to differ by “stacking flips”, which are exactly the same reshufflings that constitute our Monte Carlo tile-flips.

We also found (from simulations) that – among (small) HBS tilings in which H tiles form a network touching tip-to-tip – the in-plane tile Hamiltonian is dominated by interactions between “H(NiNi)” tiles (Hexagon tiles decorated by NiNi). Thus we suggest

$$_{\mathrm{tile}}=c_{72}N_{72}+E_{\mathrm{sf}}N_{\mathrm{sf}}.$$
(1)
where $`N_{72}`$ is the number of H(NiNi) tile pairs related by a $`72^{}`$ rotation about their common tip, and $`N_{\mathrm{sf}}`$ is the total number of stacking flips.

We fitted $`c_{72}=0.218`$eV  from the simulation with ideal sites; this is positive since the $`72^{}`$ relation of H(NiNi) tiles creates intertile NiNi neighbors, and thus reduces the number of (favorable) AlNi pairs. The structure minimizes $`N_{72}`$ by arranging that in some places, an H(NiNi) tile has one tip which contacts no other H(NiNi), and sees CoAl (at a $`144^{}`$ angle) in place of NiNi. That may be accomplished in two ways. If we re-allow lattice-gas hopping on the interior sites of the (small) HBS tiles, and properly adjust the stoichiometry , we obtain the same old rigid decoration, except that a few H(CoAl) tiles appear, as illustrated in Fig. 1. On the other hand, if we retain the rigid NiNi decoration of small H tiles, the CoAl belongs to the interior of a Boat tile, and the large HBS tiling acquires a fourth kind of tile (“Bowtie”), as in Fig. 4(a) of Ref. .

Existing experiments can in principle measure $`E_{\mathrm{sf}}`$. The time-dependent local “phason” flips , observed at $`T1200`$K by TEM through very thin samples of $`d`$(AlCuCo), can only occur by nucleating a “stacking flip” at one surface, which subsequently random-walks to the other surface. Video images  display an intermediate state for a noticeable fraction – say 10% – of the time, so that with a sample thickness estimated at $`100`$ Å (25 bilayers), one would very roughly estimate $`E_{\mathrm{sf}}0.3`$ eV in $`d`$(AlCuCo). It would be exciting if such observations could be analyzed quantitatively. In our simulation, with rigid atoms decorating ideal sites, $`E_{\mathrm{sf}}1.4`$ eV. However, relaxation (as in Fig. 2) nearly cancels this energy cost, indeed $`E_{\mathrm{sf}}0.2`$ eV in very preliminary results, which should not be compared with the experiments on a different quasicrystal. (Physically, $`E_{\mathrm{sf}}<0`$ probably implies a period $`2c8`$ Åmodulation.)

## VII Discussion

We find it remarkable that, with quite sketchy experimental input, our simulations appear competitive with single-crystal diffraction as a way of discovering this rather complex structure. This occurs despite the shortcomings of the pair potentials – our omission of 3-body terms, and the likelihood that the $`V_{\mathrm{AlTM}}`$ we use is stronger than the real one. Even if the connection to the microscopics is not quantitative, one may still obtain the large-scale order quite well: that depends mostly on the tile Hamiltonian (1) having the proper form, which is more robust than the numerical values of the coefficients.

It should also be pointed out that even imperfect potentials could be quite useful for augmenting a fit of diffraction data. A combination of energy and diffraction data could overcome minor spurious features in refinements .

Our results strongly suggest that, in this alloy system, stability of the quasicrystal requires a ternary not only to tune the electron density, as in the Hume-Rothery picture, but also because each species fills a particular type of site. Thus it would be highly desirable to explore some other composition regimes in simulations. The Co-rich “basic Co” subphase of $`d`$(AlNiCo) would be of particular interest. It exhibits much stronger diffuse scattering than Basic Ni, particularly halfway between the layers of Bragg peaks in reciprocal space, indicating a local tendency to doubling the $`c`$ periodicity ($`8`$ Å structures).

Our constraint that the atoms sit on ideal sites in the initial stage of energy optimization has considerable potential to distort the conclusions. Clearly, atom displacements have a very strong effect so it is crucial to include these for a final answer. Small adjustments of position can gain as much energy as swaps of species, and in some cases the optimal sites may lie halfway between ideal sites, or between ideal layers – the so-called “puckered” layers We do allow such relaxations in later stages, after determining a decoration, but it appears possible that a different decoration would be obtained, if we allowed relaxations in the early stages. This seems to be the most serious drawback of the lattice-gas approach  we used. Its advantages are that the systematic Monte Carlo exploration is not biased by its users’ prejudices, and furthermore it could readily be adapted to $`T>0`$ (the temperatures at which the quasicrystal phase is stable). Nonzero temperature probably affects the result just as strongly as relaxed positions do.

Within the framework of pair potentials and tile-decoration representation of the structure, a different approach is possible , whereby relaxed sites are used in the main discovery process. To implement this so that only valid tile decorations get explored, one performs relaxations with the decoration-equivalent atoms of a given class (“orbit”) being constrained to move together. However, that method is ad-hoc, depending on its users’ educated guesses as to the tilings and decorations to be tried. Future work should strive to blend the best features of these two approaches to using decorations with potentials for the discovery of structures.

This research was supported by the Department of Energy grant DE-FG02-89ER-45405 and by National Science Foundation grant DMR-0111198. We thank E. Cockayne, F. Gähler, E. Abe, and K. Edagawa for discussions.


text:
# Locality and reality revisited

## I A Lament

Many important concepts in quantum foundations are misleadingly named:

* The measurement problem is not the problem of understanding precisely what a measurement is — something which needn’t be, and probably isn’t, specified in the fundamental theory of nature — but the problem of finding a precise mathematical or physical characterisation of the sample space on which quantum probabilities are defined.
* Hidden variables are hidden in the sense that the quantum formalism conceals them, but, as Bell emphasized, certainly should not be supposed to be hidden from us, since they are meant to define the reality we experience. And Bell’s own preferred term, beable, has a metaphysically confusing connotation of a not necessarily actualised potential for existence. This would seem to make sense for the values of suitable variables — only one value is realised out of many possibilities — but not for the variables themselves: nonetheless, Bell intended the latter. Better terms for quantities parametrising the sample space which solves the measurement problem would be existents (Shimony’s term) or perhaps real variables.
* Bell’s free will criterion conflicts with some widely held views of what it means to be able to incorporate free will into a physical theory — and Bell’s response to critics explains that he did not actually intend to express a view on the problem of free will. Free will, as Bell used it, is no more than a convenient shorthand term. A theory satisfies the free will criterion if it allows us to justify treating the measurement settings as independent variables to those characterising the state of the microscopic system (typically, a pair of separated entangled particles) being measured. This need not mean that the theory admits your, or anybody’s, favourite notion of free will. Conversely, as Price emphasizes, a theory which fails to satisfy the free will criterion may nonetheless be compatible with sensible definitions of free will.
* Finally, quantum non-locality is a confusingly oxymoronic phrase. Quantum field theories — the best version of quantum theory we have — do satisfy a precise locality condition: operators with support in space-like separated regions commute. What the term is (or should be) meant to signify is that quantum theory is not locally causal under Bell’s definitions — but quantum non-(local causality) is too clumsy.

What to do? Using non-standard terms, or inventing new ones, makes life difficult for readers. So I will use the standard terms, with these criticisms understood.

## II Some motivation

It seems to me that all our difficulties in quantum foundations stem first from the measurement problem, and then from the fact that the obvious ways of tackling the measurement problem are hard to reconcile with special relativity and/or reality.

Probably this is because quantum theory is wrong, and the insolubility of the measurement problem is a symptom. Just possibly, another interpretation of quantum theory might do the trick.

Either way, there’s a case for carefully re-examining the conventional lines of argument, looking for gaps and trying to find scope for something new. Maybe there’s also a case for giving up and waiting for experiment to give us a clue, but I’ll resist that temptation here. Instead, I want to describe a different way of thinking about hidden variables, which at least sheds a different light on Bell correlations.

## III Lattice models in statistical physics

To statistical physicists, lattice models are mathematical idealizations. Many quantities of interest in statistical physics are relatively insensitive to the details of the microphysics, and many distinct systems and models may belong to the same universality class, with essentially the same properties. So, if a well chosen lattice model reproduces the bulk properties and correlation functions observed in an experimental sample, this need not imply that the model even approximately describes the system’s microphysics.

I want here to use some standard definitions of lattice statistical physics in an unfamiliar way. First, I note that lattice-like models make sense as fundamental theories in their own right, rather than as phenomenological models attempting to capture some of the properties of a more fundamental theory. Then I want to consider the possibility that something rather like a lattice model could be a fundamental theory of microphysics, giving a hidden variable theory explanation of Bell correlations.

Consider the two-dimensional Ising model on a square lattice (Figure 1). A “spin” variable $`s_i\{1,1\}`$ is attached to each site $`i`$. The Hamiltonian

$$H=\underset{(i,j)\mathrm{nearest}\mathrm{neighbours}}{}(1s_is_j),$$

includes only interactions between horizontally or vertically adjacent spins. For any possible configuration of spins $`S=\{s_i:i\mathrm{a}\mathrm{lattice}\mathrm{site}\}`$, we have a probability weight

$$p(S)=\mathrm{exp}(CH),$$

where $`C`$ is a constant (at fixed temperature). The idealised system represented by the Ising model will always be found in some legitimate spin configuration, and the probability of finding it in the configuration $`S`$ is

$$\mathrm{Prob}(S)=\frac{p(S)}{_{\mathrm{all}\mathrm{configurations}S^{}}p(S^{})}.$$

Note that the model is thus statistically local, in the sense that if two spin configurations $`S`$ and $`S^{}`$ differ only on a region $`R`$, the ratio

$$\frac{\mathrm{Prob}(S)}{\mathrm{Prob}(S^{})}$$

can be calculated purely from the values of the $`s_i`$ and $`s_i^{}`$ within the region $`R`$, together with those of their nearest neighbours. That is, if $`S`$ and $`S^{}`$ differ only locally, the ratio of their probabilities can be calculated from local quantities.

Now, the main point. The Ising model, as just described, could be interpreted as a fundamental theory. Imagine a world in which, say, a particular type of crystal has the property that its surfaces are perfect square arrays of atoms, in which each atom has an asymmetric feature which define an orientation, and in which the orientation is always perpendicular to the surface, either pointing outward ($`+1`$) or inward ($`1`$). Imagine further that the orientations are stable — we cannot change them — and can be examined by powerful microscopes, and that a plentiful supply of crystals is available.

Then it makes sense, as a hypothesis, to suppose that the atomic orientations are randomly configured, with probabilities defined by the Ising model (on the finite lattice defined by the crystal surface). The theory can be tested: examining large numbers of crystals will reveal whether or not the data are statistically consistent with the theory. And it need not be a phenomenological theory; it could be fundamental. It could simply be that nature has fixed these unchanging orientation configurations, not through some underlying dynamics whose end result is described by the Ising model, but from the beginning, by a random choice which (like the randomness of quantum measurements, according to the orthodox view) has no underlying explanation.

Obviously, this imaginary Ising-world is not much like ours, in which crystals form and decompose, atoms don’t behave classically, their configurations and spin states aren’t fixed, and everything is (as far as we presently understand) ultimately described by quantum dynamics. But conceptually the Ising theory of crystals in Ising-world makes sense. It is testable in much the same way as quantum theory is, and it includes a simplifying principle — statistical locality — which makes it much easier to test: it suffices to look at orientation configurations in local regions of crystal surfaces. And, given our other assumptions about Ising-world, there is no compelling reason — nothing analogous to the measurement problem — to think Ising theory must be incomplete.

## IV Minkowski space trajectory models

Let us suppose that underlying quantum theory there is some theory in which real particles follow definite trajectories through spacetime. To simplify things, take space-time to be two-dimensional, with a discrete lattice structure at microscopic scales (Figure 2). The lattice edges are lightlike, and so the discussion is restricted to massless particles, which travel along the lattice edges. Also, for the purposes of this discussion, macroscopic experimental apparatus will be treated as qualitatively distinct from microscopic particles, so that we are allowed to consider apparatus set up somewhere in spacetime without having to describe it as a collection of microscopic particles, each following their own trajectories through the lattice.

All of these simplifications are meant to be provisional rather than fundamental. The idea is to set up the standard discussion of Bell correlations, and then find a new way of looking at them, without obscuring things with too many technicalities.

Figure 3 shows a standard Bell experiment on a pair of entangled massless particles, “photons”, with one extra feature — lightlike trajectories along the lattice edges have been included. Though this looks like a standard beam path picture, the trajectories here are hidden variables: the picture is meant to be interpreted as saying that the photons really did travel along the paths depicted — something we will eventually have to show is consistent with the model we produce. The space-time trajectories of the classical measuring apparatuses won’t affect the discussion, and are omitted. To further simplify things — since we cannot use photon polarizations in two dimensions — let us assume that the measurements are perfect von Neumann measurements of an internal degree of freedom, “spin”, with the photons initially in the singlet state of a pair of spin $`1/2`$ particles.***The scare quotes around “photon” and “spin” are left implicit from now on. According to quantum theory, these von Neumann measurements do not affect the form of the photon wave packet, assuming (let’s) that the spatial and spin degrees of freedom are unentangled. It seems reasonable to postulate that, in the hidden variable model we are building, the photon trajectories will also be unaffected by such measurements. For instance, after a measurement of spin in the measurement basis defined by vector $`𝐚`$, the apparatus records the answer $`𝐚`$ or $`𝐚`$, the photon emerges in the corresponding spin eigenstate, and its trajectory continues unaltered.

The aim now is to reproduce Bell correlations by a statistically local lattice-like model involving photon trajectories, while respecting Bell’s free will criterion. We will want to attach further hidden variable degrees of freedom to the photon trajectories, so that the spin measurements become relevant to the hidden variable model. The free will criterion requires that no hidden variables at points before the measurements should be correlated with the measurement choices. This, together with locality, means that the photons’ states before the measurements cannot explain the correlations: we need rules involving events after the measurements.

We can ensure a relevant post-measurement event if the photons are forced to meet again, and that can be done, without violating statistical locality, as follows. First, we give a small but non-zero probability weight to trajectory segments in which a photon’s direction switches: this allows the two photons to meet after the measurements. Second, we choose rules that assign probability zero to the event that the photons propagate forever, and non-zero probability weight to the event that two photon trajectories annihilate when they meet. Together these ensure that photons must eventually undergo pair annihilation. The first of these can be arranged by assigning probability weight of $`(1ϵ)`$ to a vertex in which the photon continues in the same direction, and weight $`ϵ`$ to a vertex in which it switches (Figure 4); here $`ϵ>0`$ is chosen to be small enough that $`(1ϵ)^N1`$, where $`N`$ is the number of lattice links a photon would traverse in the largest scale Bell experiment carried out to date.

In choosing rules for a pair annihilation, we need to make sure that the originally entangled photons meet and mutually annihilate each other, and not other photons. This we do by introducing a hidden variable label which only the entangled pair share, and allowing annihilation only between pairs which share the same label. Finally, we enforce the correct correlations by fiat, by adding a further label which records a photon’s measurement outcome, and by making the probability weight for the annihilation event depend appropriately on these labels (Figure 5). To summarise: we’ve arranged rules so that the entangled photons from pair $`i`$ do the following: (i) with high probability, propagate linearly from source to detectors, (ii) after measurement, each eventually undergo one or more changes of direction, (iii) with probability one, eventually meet up and annihilate. (Figure 6)

## V Discussion

All we’ve done, of course, is rewrite the standard quantum correlations into ad hoc rules, in such a way that we can formally attach hidden variables characterising the measurement results to the photon paths after measurement, and formally account for the correlations as deriving from a rule attached to a pair annihilation event that must take place in the far future. What was the point of that?

The point, as I see it, is to underline how difficult it is to draw absolutely watertight conclusions about the form of realist theories underlying quantum theory, even if we accept that such theories should respect special relativity and should incorporate some form of locality, and even after decades of analysis of the implications of Bell experiments. To recapitulate:

We can agree with Bell that any theory which postulates common causes in the past which systematically affect experimenter and system, in such a way as to always produce Bell correlations, is likely to be unbelievably conspiratorial. But, as Price has so carefully and persuasively argued, this leaves open the possibility of theories which invoke some form of advanced action principle, by which experimenters’ choices affect the earlier state of the system. As Price notes, although such theories would violate Bell’s free will criterion, they would not necessarily rule out free will per se, at least according to some reasonable definitions. And they would not contradict the notion of effective freedom which underlay Bell’s definition, nor invoke the same type of implausible conspiracy as a common cause explanation.

Still, there is (at the very least) an aesthetic case for the free will criterion. But even if the criterion is accepted, the discussion isn’t closed. The models described above respect the free will criterion and they invoke no complex common cause hypothesis. They incorporate a respectable form of locality. Yet they give — in the simple case considered — a consistent explanation of Bell experiments.

That said, there is, of course, much to be said against the sort of ideas outlined above. The rules proposed are entirely ad hoc, transparently constructed precisely in order to reproduce quantum correlations. They could be generalised to cover more complicated entanglement experiments, but only — it would seem — at the price of introducing even more ad hoc and unwieldy rules. They require, in particular, that every hidden particle trajectory carries a record of its past history — which particles it was originally entangled with, and which measurements were carried out on it — arbitrarily far into the future. It is hard to see, methodologically speaking, what the point of pursuing this idea would be, unless some independent principle could be found to explain why the annihilation vertex weights should take the form necessary to reproduce quantum correlations. And where could such a principle come from?

In fact, any type of correlation, including those which allow superluminal signalling, could consistently be attached to the annihilation vertex. Why should such correlations be excluded? If we were to take the basic idea of these models seriously, we would be led to ask why the non-local correlations we see in nature are so tame as to preclude superluminal signalling. That this would lead to paradoxes is not a valid answer. Consistent probabilistic rules for global configurations cannot produce a logical contradiction. Attempts to set up a paradoxical causal loop fail. This is essentially because the probability of a successful superluminal signal is — even in isolation — always at least slightly less than one, and the probability of a sequence of signals succeeding is not generally the product of the values of the individual probabilities of success if the relevant signalling devices were isolated.

Moreover, the model involved several breathtaking simplifications, which would need to be removed if the idea were to be taken seriously. (Can we eliminate the need for an underlying lattice, or ensure that the model has a Lorentz invariant continuum limit? What about rules for massive particles, and for other types of measurement? And how do we build a coherent model that describes macroscopic measuring devices in terms of microscopic particles?)

Nonetheless, despite all these flaws, I think there is some interest in this rewriting of the standard story. It seems to show that the conventional wisdom about Bell’s “free will” criterion is not quite right: realist, local and Lorentz invariant theories can respect the criterion and yet reproduce quantum correlations. And it does so by identifying a class of realist models — probabilistic models in which the probability of a configuration depends on a product of local terms — which standard discussions of Bell’s theorem neglect. There are, as just noted, many arguments against taking seriously the idea that nature is actually described by these models. But the arguments are different from the arguments against other unorthodox ideas about hidden variables and Bell correlations, and worth thinking through, just for that reason.

## VI Acknowledgements

I am very grateful to Tomasz Placek and Jeremy Butterfield for giving me the opportunity to present this work at a most stimulating and enjoyable meeting, and for patient editorial encouragement far beyond the call of duty. The work was supported in part by a Royal Society University Research Fellowship, PPARC, and the European collaboration EQUIP.


text:
# Variations in the Bivariate Brightness Distribution with different galaxy types.

## 1. Introduction

The BBD is the number density of galaxies as a function of both absolute magnitude, $`M`$, and effective surface brightness, $`\mu _e`$. This provides a useful framework for removing surface brightness selection effects from the measurement of the space density of galaxies (Cross et al. 2001) and also provides useful constraints for testing models of galaxy formation and evolution (de Jong & Lacey 2000).

Large redshift surveys such as the “Two-degree Field Galaxy Redshift Survey” (2dFGRS, Colless et al. 2001) and the “Sloan Digital Sky Survey (SDSS, York et al. 2000) provide the quantity of data necessary to select the unbiased, statistically significant samples needed to measure the BBDs of different galaxy types.

## 2. The Data

The 2dFGRS is a $`2000_{}^{2}{}_{}{}^{}`$ redshift survey, containing redshifts for $`220,000`$ galaxies, with an apparent magnitude limit $`b_j=19.45`$ (see Colless et al. 2001). The input catalogue for the 2dFGRS is based on the Automated Plate Machine (APM) catalogue (Maddox et al. 1990), with some recent recalibration using new CCD data.

Madgwick et al. (2002) used the 2dFGRS to identify four spectral classes based on the emission/absorption line strength ($`\eta `$) of galaxies. $`\eta `$-type 1 galaxies have strong absorption lines and $`\eta `$-type 4 galaxies have strong emission lines.

The Millennium Galaxy Catalogue (MGC, Liske et al. 2002) is a $`36`$deg2 deep CCD imaging survey of an equatorial strip between $`10^h`$ and $`15^h`$. The MGC images were taken through the KPNO-B filter on the Isaac Newton Telescope. The detection isophote $`\mu _{lim}=26.0`$ mag arcsec-2 and the mean seeing is $`1.3^{\prime \prime }`$. Star-galaxy separation has been undertaken by eye for all $`B_{MGC}<20`$ objects. For each galaxy we measured half light radii ($`r_e`$) along the major axis of the ellipse that contains half of the flux of the galaxy. We use this to calculate the effective surface brightness, $`\mu _e=B_{MGC}+2.5\mathrm{log}_{10}(2\pi r_e^2)`$, assuming a circular aperture, to correct for inclination effects, see Cross et al. (2002).

## 3. Photometry and Completeness of the 2dFGRS

Using a mean colour for galaxies $`(\overline{BV})=1.0`$, we find that the filter conversion between the APM and MGC is $`B_{MGC}b_j=0.1353(BV)=0.1353`$. $`\mathrm{\Delta }m(MGC2dFGRS)=0.078\pm 0.157`$.

The mean scale error in magnitude is small and linear (see Fig 1). It is fit by the equation $`\mathrm{\Delta }m=a+b(B_{MGC}19.45)`$, where $`a=0.0467`$ and $`b=0.0221`$. However, there is a significant scale error as a function of effective surface brightness (see Fig 1). $`\mathrm{\Delta }m`$ is constant for $`22<\mu _e<25`$ but decreases non-linearly with $`\mu _e`$ for $`\mu _e<22`$ mag arcsec-2. The most likely explanation is non-linearities or saturation of high surface brightness galaxies in the APM. A comparison between the MGC and the Sloan Digital Sky Survey Early Data Release (SDSS-EDR) shows that $`\mathrm{\Delta }m(MGCSDSS)=0.009\pm 0.100`$ with no significant scale error.

We select the sample of MGC galaxies which should have 2dFGRS matches to measure the photometric completeness, see Cross et al. (2002). The variation of incompleteness with magnitude and with surface brightness is shown in Fig. 1. The incompleteness is worst at the bright and faint ends, reaching a minimum of $`5\%`$ in the range $`17.25<B_{MGC}<18.0`$. It increases linearly with magnitude in the range $`18.0<B_{MGC}<19.0`$. For $`B_{MGC}>19.0`$, the magnitude limit of the 2dFGRS affects the measurement of completeness. For $`22.25<\mu _e<24.25`$ the incompleteness is fairly constant $`5\%`$. The incompleteness of low surface brightness galaxies increases rapidly beyond $`\mu _e=24.25`$ and no 2dFGRS galaxies are seen with $`\mu _e>25.75`$, as expected with an isophotal limit $`\mu _{lim}=24.67`$. At the high surface brightness end, the incompleteness rises steadily, since small angular size galaxies were mistaken for stars in the APM.

## 4. Bivariate Brightness Distributions

Since the 2dFGRS suffers from a non-linearity with surface brightness and $`8\%`$ photometric incompleteness to $`B_{MGC}>19.0`$, we use the MGC objects defined as galaxies to select the BBD sample. The selection limits used are: the isophotal limit $`\mu _{lim,B}=26.0`$ mag arcsec-2, the minimum half light radius $`r_e=0.793^{\prime \prime }`$, defined by the average star-galaxy separation, and the maximum half-light radius $`r_e=15.8^{\prime \prime }`$, limited by the sky subtraction. The other limits, $`15.0<B_{MGC}<19.0,21.0<\mu _e<24.5`$ are determined by the redshift completeness. We select a subsample, with $`>70\%`$ redshift completeness at all values of $`B_{MGC}`$ and $`\mu _e`$. These limits are discussed in Cross (2002).

We calculate the overall BBD using a maximum likelihood method outlined in Cross (2002). It is based on the step-wise method (Efstathiou, Ellis & Peterson 1988), but also includes Visibility Theory (Phillipps, Davies & Disney 1990). The overall BBD is shown in the first plot of Fig 2.

$$\varphi (M,\mu _e)=\frac{0.4\mathrm{ln}(10)}{\sqrt{2\pi }\sigma _\mu }\varphi _{}\mathrm{\hspace{0.17em}10}^{0.4(M^{}M)(\alpha +1)}e^{10^{0.4(M^{}M)}}e^{\left[\frac{1}{2}\left(\frac{\mu _e\mu _e^{}\beta _\mu (MM_{})}{\sigma _\mu }\right)^2\right]}$$
(1)
We fit the Chołoniewski (1985) function, Eqn 1, to this distribution. The best fit gives $`M^{}=19.24`$, $`\varphi ^{}=2.47\times 10^2`$, $`\alpha =0.796`$, $`\beta _\mu =0.180`$, $`\mu _e^{}=22.06`$ and $`\sigma _\mu =0.765`$ with $`\chi ^2=220`$ for $`\nu =101`$ degrees of freedom.

The MGC database is too small to produce statistically significant samples to measure the BBD for each spectral type. Instead, we use the MGC to remove the non-linearities in the 2dFGRS magnitudes and to correct for photometric incompleteness and stellar contamination, see Cross (2002). This gives a ten times larger dataset, but with larger random errors in the photometry.

The BBDs produced for each type are also displayed in Fig 2. We find that $`\eta `$-type 1 galaxies have a bimodal distribution - a tightly bounded distribution, centred on $`M=18.5,\mu _e=22.2`$ with a high degree of symmetry and no significant luminosity surface brightness correlation, with a second fainter distribution, which continues upto the selection limits.

In contrast types 2-4 show a strong luminosity-surface brightness correlation ($`\beta _\mu =0.25\pm 0.05`$, $`\sigma _\mu =0.56\pm 0.01`$), which does not vary significantly between these types. The proportion of faint galaxies increases with $`\eta `$-type: $`\eta `$-type 2 galaxies have a flat faint end slope ($`\alpha =0.8`$) and $`\eta `$-type 4 galaxies have a steep faint end slope ($`\alpha =1.5`$).

## 5. Conclusions

We have used a deep, wide field CCD survey, the MGC, to show that the 2dFGRS contains a scale error with surface-brightness, due to non-linearities in the photographic plates. Furthermore, the 2dFGRS has a mean photometric incompleteness of $`8\%`$, for $`B_{MGC}<19`$, with many low surface brightness and low angular resolution galaxies excluded.

The MGC can be used to define a dataset with well understood selection limits that contains CCD photometric data and 2dFGRS spectroscopic data. We have used this dataset to calculate the bivariate brightness distribution (BBD) for the overall galaxy distribution and for four spectral types.

We show that strongly absorbing galaxies (ellipticals/lenticulars) form a bimodal population. The bright, high surface brightness population does not show a strong correlation between magnitude and surface brightness. $`\eta `$-type 2-4 galaxies (spirals/irregulars) all show the same strong correlation between magnitude and surface brightness. The proportion of dwarf galaxies increases with increasing emission line strength.

## References

Chołoniewski J. 1985, MNRAS, 214, 197

Colless M. et al. 2001, MNRAS, in press

Cross N.J.G. et al. 2001, MNRAS, 324, 825

Cross N.J.G., Driver S.P. 2002, MNRAS, 329, 579

Cross N.J.G. et al. 2002, in preparation

Cross N.J.G. 2002, PhD Thesis, Univ. St Andrews.

de Jong R., Lacey C. 2000, ApJ, 545, 781

Efstathiou G., Ellis R., Peterson B. 1988, MNRAS, 232, 431

Liske J., Lemon D.J., Driver S.P., Cross N.J.G., Couch W.J.C. 2002, in preparation

Maddox S.J., Sutherland W.J., Efstathiou G., Loveday J. 1990, MNRAS, 246, 433

Madgwick D. et al. 2002, MNRAS, accepted

Phillips S., Davies J., Disney M. 1990, MNRAS, 242, 235

York D.G. et al. 2000, AJ, 120, 1579