text:
# A Combinatorial Interpretation of The Numbers 6(2𝑛)!/𝑛!(𝑛+2)!

## 1. Introduction

The Catalan numbers

$$C_n=\frac{1}{n+1}\left(\genfrac{}{}{0pt}{}{2n}{n}\right)=\frac{(2n)!}{n!(n+1)!}$$

are well-known integers that arise in many combinatorial problems. Stanley \[9, pp. 219–229\], gives $`66`$ combinatorial interpretations of these numbers.

In 1874 E. Catalan  observed that the numbers

(1.1)
$$\frac{(2m)!(2n)!}{m!n!(m+n)!}$$

are integers, and their number-theoretic properties were studied by several authors (see Dickson \[2, pp. 265–266\]). For $`m=0`$, (1.1) is the middle binomial coefficient $`\left(\genfrac{}{}{0pt}{}{2n}{n}\right)`$, and for $`m=1`$ it is $`2C_n`$.

Except for $`m=n=0`$, these integers are even, and it is convenient for our purposes to divide them by 2, so we consider the numbers

(1.2)
$$T(m,n)=\frac{1}{2}\frac{(2m)!(2n)!}{m!n!(m+n)!}.$$

Some properties of these numbers are given in , where they are called “super Catalan numbers”. An intriguing problem is to find a combinatorial interpretation to the super Catalan numbers. The following identity \[3, Equation (32)\], together with the symmetry property $`T(m,n)=T(n,m)`$ and the initial value $`T(0,0)=1`$, shows that $`T(m,n)`$ is a positive integer for all $`m`$ and $`n`$.

(1.3) $`{\displaystyle \underset{n}{}}2^{p2n}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{2n}}\right)T(m,n)=T(m,m+p),p0.`$
Formula (1.3) allows us to construct recursively a set of cardinality $`T(m,n)`$ but we have not found any natural description of it for $`m2`$. Shapiro  gave a combinatorial interpretation to (1.3) in the case $`m=1`$, which is the Catalan number identity

$$\underset{n}{}2^{p2n}\left(\genfrac{}{}{0pt}{}{p}{2n}\right)C_n=C_{p+1}.$$

A similar interpretation works for the case $`m=0`$ of (1.3) (when multiplied by 2), which is the identity

$$\underset{n}{}2^{p2n}\left(\genfrac{}{}{0pt}{}{p}{2n}\right)\left(\genfrac{}{}{0pt}{}{2n}{n}\right)=\left(\genfrac{}{}{0pt}{}{2p}{p}\right).$$

Another intriguing formula for $`T(m,n)`$, which does not appear in , is

(1.4)
$$1+\underset{m,n=1}{\overset{\mathrm{}}{}}C_mC_nx^my^n=\left(1\underset{m,n=1}{\overset{\mathrm{}}{}}T(m,n)x^my^n\right)^1.$$

Although (1.4) suggests a combinatorial interpretation for $`T(m,n)`$ based on a decomposition of pairs of objects counted by Catalan numbers, we have not found such an interpretation.

In this paper, we give a combinatorial interpretation for $`T(2,n)=6(2n)!/n!(n+2)!`$ for $`n1`$ and for $`T(3,n)=60(2n)!/n!(n+3)!`$ for $`n2`$. The first few values of $`T(m,n)`$ for $`m=2`$ and $`m=3`$ are as follows:

| $`m\backslash n`$ | $`0`$ | $`1`$ | $`2`$ | $`3`$ | $`4`$ | $`5`$ | $`6`$ | $`7`$ | $`8`$ | $`9`$ | $`10`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`2`$ | $`3`$ | $`2`$ | $`3`$ | $`6`$ | $`14`$ | $`36`$ | $`99`$ | $`286`$ | $`858`$ | $`2652`$ | $`8398`$ |
| $`3`$ | $`10`$ | $`5`$ | $`6`$ | $`10`$ | $`20`$ | $`45`$ | $`110`$ | $`286`$ | $`780`$ | $`2210`$ | $`6460`$ |

We show that $`T(2,n)`$ counts pairs of Dyck paths of total length $`2n`$ with heights differing by at most 1. We give two proofs of this result, one combinatorial and one using generating functions. The combinatorial proof is based on the easily checked formula

(1.5)
$$T(2,n)=4C_nC_{n+1}$$

which we interpret by inclusion-exclusion.

Our interpretation for $`T(3,n)`$ is more complicated, and involves pairs of Dyck paths with height restrictions. Although we have the formula $`T(3,n)=16C_n8C_{n+1}+C_{n+2}`$ analogous to $`(\text{1.5})`$, we have not found a combinatorial interpretation to this formula, and our proof uses generating functions.

Interpretations of the number $`T(2,n)`$ in terms of trees, related to each other, but not, apparently, to our interpretation, have been found by Schaeffer , and by Pippenger and Schleich \[5, pp. 34\].

## 2. The Main theorem

All paths in this paper have steps $`(1,1)`$ and $`(1,1)`$, which we call *up steps* and *down steps*. A step from a point $`u`$ to a point $`v`$ is denoted by $`uv`$. The *level* of a point in a path is defined to be its $`y`$-coordinate. A *Dyck path* of *semilength* $`n`$ (or of length $`2n`$) is a path that starts at $`(0,0)`$, ends at $`(2n,0)`$, and never goes below level $`0`$. It is well-known that the number of Dyck paths of semilength $`n`$ equals the Catalan number $`C_n`$. The *height* of a path $`P`$, denoted by $`h(P)`$, is the highest level it reaches.

Every nonempty Dyck path $`R`$ can be factored uniquely as $`UPDQ`$, where $`U`$ is an up step, $`D`$ is a down step, and $`P`$ and $`Q`$ are Dyck paths. Thus the map $`P(P,Q)`$ is a bijection from nonempty Dyck paths to pairs of Dyck paths. Let $`𝐁_n`$ be the set of pairs of Dyck paths $`(P,Q)`$ of total semilength $`n`$. This bijection gives $`|𝐁_n|=C_{n+1}`$, so by (1.5), we have $`T(2,n)=4C_n|𝐁_n|`$.

Our interpretatoin for $`T(2,n)`$ is a consequence of the following Lemma 2.1. We give two proofs of this lemma, one combinatorial and the other algebraic. The algebraic proof will be given in the next section.

###### Lemma 2.1.

For $`n1`$, $`C_n`$ equals the number of pairs of Dyck paths $`(P,Q)`$ of total semilength $`n`$, with $`P`$ nonempty and $`h(P)h(Q)+1`$.

###### Proof.

Let $`𝐃_n`$ be the set of Dyck paths of semilength $`n`$, and let $`𝐄_n`$ be the set of pairs of Dyck paths $`(P,Q)`$ of total semilength $`n`$, with $`P`$ nonempty and $`h(P)h(Q)+1`$.

For a given pair $`(P,Q)`$ in $`𝐄_n`$, since $`P`$ is nonempty, the last step of $`P`$ must be a down step, say, $`uv`$. By replacing $`uv`$ in $`P`$ with an up step $`uv^{}`$, we get a path $`F_1`$. Now raising $`Q`$ by two levels, we get a path $`F_2`$. Thus $`F:=F_1F_2`$ is a path that ends at level $`2`$ and never goes below level $`0`$. The point $`v^{}`$ belongs to both $`F_1`$ and $`F_2`$, but we treat it as a point only in $`F_2`$, even if $`F_2`$ is the empty path. The condition that $`h(P)h(Q)+1`$ yields $`h(F_1)<h(F_2)`$, which implies that the highest point of $`F`$ must belong to $`F_2`$. See Figure 1 below.

Now let $`y`$ be the leftmost highest point of $`F`$ (which is in $`F_2`$), and let $`xy`$ be the step in $`F`$ leading to $`y`$. Then $`xy`$ is an up step. By replacing $`xy`$ with a down step $`xy^{}`$, and lowering the part of $`F_2`$ after $`y`$ by two levels, we get a Dyck path $`D𝐃_n`$. See Figure 2 below.

With the following two key observations, it is easy to see that the above procedure gives a bijection from $`𝐄_n`$ to $`𝐃_n`$. First, $`x`$ in the final Dyck path $`D`$ is the rightmost highest point. Second, $`u`$ in the intermediate path $`F`$ is the rightmost point of level $`1`$ in both $`F`$ and $`F_1`$. ∎

###### Theorem 2.2.

For $`n1`$, the number $`T(2,n)`$ counts pairs of Dyck paths $`(P,Q)`$ of total semilength $`n`$ with $`|h(P)h(Q)|1`$.

###### Proof.

Let $`F`$ be the set of pairs of Dyck paths $`(P,Q)`$ with $`h(P)h(Q)+1`$, and let $`G`$ be the set of pairs of Dyck paths $`(P,Q)`$ with $`h(Q)h(P)+1`$. By symmetry, we see that $`|F|=|G|`$. Now we claim that the cardinality of $`F`$ is $`2C_n`$. This claim follows from Lemma 2.1 and the fact that if $`P`$ is the empty path, then $`h(P)h(Q)+1`$ for every $`Q𝐃_n`$.

Clearly we have that $`FG=𝐁_n`$, and that $`FG`$ is the set of pairs of Dyck paths $`(P,Q)`$, with $`|h(P)h(Q)|1`$. The theorem then follows from the following computation:

$`|FG|=|F|+|G||FG|=4C_n|𝐁_n|=4C_nC_{n+1}.`$
## 3. An Algebraic Proof and Further Results

In this section we give an algebraic proof of Lemma 2.1.

Let $`c(x)`$ be the generating function for the Catalan numbers, so that

$$c(x)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{n+1}\left(\genfrac{}{}{0pt}{}{2n}{n}\right)x^n=\frac{1\sqrt{14x}}{2x}.$$

Then $`c(x)`$ satisfies the functional equation $`c(x)=1+xc(x)^2`$. Let $`C=xc(x)^2=c(x)1`$ and let $`G_k`$ be the generating function for Dyck paths of height at most $`k`$. Although $`G_k`$ is a rational function, a formula for $`G_k`$ in terms of $`C`$ will be of more use to us than the explicit formula for $`G_k`$.

###### Lemma 3.1.

For $`k1`$,

(3.1)
$$G_k=(1+C)\frac{1C^{k+1}}{1C^{k+2}}.$$

###### Proof.

Let $`P`$ be a path of height at most $`k1`$. If $`P`$ is nonempty then $`P`$ can be factored as $`UP_1DP_2`$, where $`U`$ is an up step, $`P_1`$ is a Dyck path of height at most $`k1`$ (shifted up one unit), $`D`$ is a down step, and $`P_2`$ is a Dyck path of height at most $`k`$. Thus $`G_k=1+xG_{k1}G_k`$, so $`G_k=1/(1xG_{k1})`$. Equation (3.1) clearly holds for $`k=1`$ and $`k=0`$. Now suppose that for some $`k1`$,

$$G_{k1}=(1+C)\frac{1C^k}{1C^{k+1}}.$$

Then the recurrence, together with the formula $`x=C/(1+C)^2`$, gives

$`G_k`$ $`=\left[1x(1+C){\displaystyle \frac{1C^k}{1C^{k+1}}}\right]^1=\left[1{\displaystyle \frac{C}{1+C}}{\displaystyle \frac{1C^k}{1C^{k+1}}}\right]^1`$
$`=\left[{\displaystyle \frac{1C^{k+2}}{(1+C)(1C^{k+1})}}\right]^1=(1+C){\displaystyle \frac{1C^{k+1}}{1C^{k+2}}}.`$
We can prove Lemma 2.1 by showing that $`_{n=0}^{\mathrm{}}G_{n+1}(G_nG_{n1})=1+2C`$; this is equivalent to the statement that the number of pairs $`(P,Q)`$ of Dyck paths of semilength $`n>0`$ with $`h(P)h(Q)+1`$ is $`2C_n`$.

###### Theorem 3.2.

$$\underset{n=0}{\overset{\mathrm{}}{}}(G_nG_{n1})G_{n+1}=1+2C.$$

###### Proof.

Let

$$\mathrm{\Psi }_k=\underset{n=k}{\overset{\mathrm{}}{}}\frac{C^n}{1C^n}.$$

Thus if $`j<k`$ then

(3.2)
$$\mathrm{\Psi }_j=\underset{n=j}{\overset{k1}{}}\frac{C^n}{1C^n}+\mathrm{\Psi }_k.$$

We have

$$G_nG_{n+1}=(1+C)^2\frac{1C^{n+1}}{1C^{n+3}}$$

and

$$G_{n1}G_{n+1}=(1+C)^2\frac{(1C^n)(1C^{n+2})}{(1C^{n+1})(1C^{n+3})}.$$

Let

$$S_1=\underset{n=0}{\overset{\mathrm{}}{}}\left(\frac{1C^{n+1}}{1C^{n+3}}1\right)$$

and

$$S_2=\underset{n=0}{\overset{\mathrm{}}{}}\left(\frac{(1C^n)(1C^{n+2})}{(1C^{n+1})(1C^{n+3})}1\right).$$

Then $`_{n=0}^{\mathrm{}}(G_nG_{n1})G_{n+1}=(1+C)^2(S_1S_2)`$. We have

$$\frac{1C^{n+1}}{1C^{n+3}}1=\frac{(1C^2)C^{n+1}}{1C^{n+3}},$$

so $`S_1=(1C^2)C^2\mathrm{\Psi }_3`$, and

$$\frac{(1C^n)(1C^{n+2})}{(1C^{n+1})(1C^{n+3})}1=\frac{(1C)C^n}{(1+C)(1C^{n+1})}\frac{(1C^3)C^{n+1}}{(1+C)(1C^{n+3})},$$

so

$$S_2=\frac{1C}{1+C}C^1\mathrm{\Psi }_1\frac{1C^3}{1+C}C^2\mathrm{\Psi }_3.$$

Therefore

$`S_1S_2`$ $`={\displaystyle \frac{1C}{1+C}}C^1\mathrm{\Psi }_1+\left({\displaystyle \frac{1C^3}{1+C}}(1C^2)\right)C^2\mathrm{\Psi }_3`$
$`={\displaystyle \frac{1C}{1+C}}C^1(\mathrm{\Psi }_1\mathrm{\Psi }_3)`$
$`={\displaystyle \frac{1C}{1+C}}C^1\left({\displaystyle \frac{C}{1C}}+{\displaystyle \frac{C^2}{1C^2}}\right)={\displaystyle \frac{1+2C}{(1+C)^2}}.`$
Thus $`(1+C)^2(S_1S_2)=1+2C`$. ∎

By similar reasoning, we could prove Theorem 2.2 directly: The generating function for pairs of paths with heights differing by at most 1 is

$$\underset{n=0}{\overset{\mathrm{}}{}}(G_nG_{n1})(G_{n+1}G_{n2}),$$

where we take $`G_1=G_2=0`$, and a calculation like that in the proof of Theorem 3.2 shows that this is equal to

$`1+2CC^2`$ $`=4c(x)c(x)^22=4c(x){\displaystyle \frac{c(x)1}{x}}2`$
$`=1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(4C_nC_{n+1})x^n=1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}T(2,n)x^n.`$
Although the fact that the series in Theorem 3.2 telescopes may seem surprising, we shall see in Theorem 3.4 that it is a special case of a very general result on sums of generating functions for Dyck paths with restricted heights.

###### Lemma 3.3.

Let $`R(z,C)`$ be a rational function of $`z`$ and $`C`$ of the form

$$\frac{zN(z,C)}{_{i=1}^m(1zC^{a_i})},$$

where $`N(z,C)`$ is a polynomial in $`z`$ of degree less than $`m`$, with coefficients that are rational functions of $`C`$, and the $`a_i`$ are distinct positive integers. Let $`L=lim_z\mathrm{}R(z,C)`$. Then

$$\underset{n=0}{\overset{\mathrm{}}{}}R(C^n,C)=Q(C)+L\mathrm{\Psi }_1,$$

where $`Q(C)`$ is a rational function of $`C`$.

###### Proof.

First we show that the lemma holds for $`R(z,C)=z/(1zC^a)`$. In this case, $`L=lim_z\mathrm{}R(z,C)=C^a`$ and

$$\underset{n=0}{\overset{\mathrm{}}{}}R(C^n,C)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{C^n}{1C^{n+a}}=C^a\underset{n=a}{\overset{\mathrm{}}{}}\frac{C^n}{1C^n}=\underset{n=0}{\overset{a1}{}}\frac{C^{na}}{1C^n}+C^a\mathrm{\Psi }_1.$$

Now we consider the general case. Since $`R(z,C)/z`$ is a proper rational function of $`z`$, it has a partial fraction expansion

$$\frac{1}{z}R(z,C)=\underset{i=1}{\overset{m}{}}\frac{U_i(C)}{1zC^{a_i}}$$

for some rational functions $`U_i(C)`$, so

$$R(z,C)=\underset{i=1}{\overset{m}{}}U_i(C)\frac{z}{1zC^{a_i}}.$$

The general theorem then follows by applying the special case to each summand. ∎

###### Theorem 3.4.

If $`i_1,i_1^{},\mathrm{},i_m`$ are distinct integers, then

$$\underset{n=0}{\overset{\mathrm{}}{}}(G_{n+i_1}G_{n+i_1^{}})G_{n+i_2}G_{n+i_3}\mathrm{}G_{n+i_m}$$

is a rational function of $`C`$.

###### Proof.

Apply Lemma 3.3 to

$$\underset{n=0}{\overset{\mathrm{}}{}}\left(G_{n+i_1}G_{n+i_2}G_{n+i_3}\mathrm{}G_{n+i_m}1\right)$$

and

$$\underset{n=0}{\overset{\mathrm{}}{}}\left(G_{n+i_1^{}}G_{n+i_2}G_{n+i_3}\mathrm{}G_{n+i_m}1\right).$$

## 4. A combinatorial interpretation for $`T(3,n)`$

It is natural to ask whether there are combinatorial interpretations to $`T(m,n)`$ for $`m>2`$ similar to Theorem 2.2. It is straightforward (with the help of a computer algebra system) to evaluate sums like (3.4) that count $`m`$-tuples of paths with height restrictions, and we find that sums like that in Theorem 3.4 involving products of $`m`$ path generating functions may generally be expressed in the form

$$R_1(x)+R_2(x)\underset{n=0}{\overset{\mathrm{}}{}}T(m,n)x^n,$$

where $`R_1(x)`$ and $`R_2(x)`$ are rational functions of $`x`$. However, for general $`m`$ we have not found a set of $`m`$-tuples of paths counted by $`T(m,n)`$. We have found a set of paths counted by $`T(3,n)`$, though it is not as simple as one would like.

We need to consider paths that end at levels greater than 0. Let us define a *ballot path* to be a path that starts at level 0 and never goes below level 0.

In the previous section all our paths had an even number of steps, so it was natural to assign a path with $`n`$ steps the weight $`x^{n/2}`$. We shall continue to weight paths in this way, even though some of our paths now have odd lengths.

Let $`G_k^{(j)}`$ be the generating function for ballot paths of height at most $`k`$ that end at level $`j`$.

###### Lemma 4.1.

For $`0jk+1`$ we have

(4.1)
$$G_k^{(j)}=C^{j/2}(1+C)\frac{1C^{kj+1}}{1C^{k+2}}$$

###### Proof.

The case $`j=0`$ is Lemma 3.1. Now let $`W`$ be a ballot path counted by $`G_k^{(j)}`$, where $`j>0`$, so that $`W`$ is of height at most $`k`$ and $`W`$ ends at level $`j`$. Then $`W`$ can be factored uniquely as $`W_1UW_2`$, where $`W_1`$ is a path of height at most $`k`$ that ends at level 0 and $`W_2`$ is a path from level 1 to level $`j`$ that never goes above level $`k`$ nor below level 1. Using $`\sqrt{x}=\sqrt{C/(1+C)^2}=\sqrt{C}/(1+C)`$, we obtain

$$G_k^{(j)}=G_k\sqrt{x}G_{k1}^{(j1)}=(1+C)\frac{1C^{k+1}}{1C^{k+2}}\frac{\sqrt{C}}{1+C}G_{k1}^{(j1)}=\sqrt{C}\frac{1C^{k+1}}{1C^{k+2}}G_{k1}^{(j1)},$$

and (4.1) follows by induction.

We note an alternative formula that avoids half-integer powers of $`C`$,

$$G_k^{(j)}=x^{j/2}(1+C)^{j+1}\frac{1C^{kj+1}}{1C^{k+2}},$$

which follow easily from (4.1) and the formula $`\sqrt{x}=\sqrt{C}/(1+C)`$.

Although we will not need it here, there is a similar formula for the generating function $`G_k^{(i,j)}`$ for paths of height at most $`k`$ that start at level $`i`$, end at level $`j`$, and never go below level 0: for $`0ijk+1`$ we have

(4.2)
$$G_k^{(i,j)}=C^{(ji)/2}(1+C)\frac{(1C^{i+1})(1C^{kj+1})}{(1C)(1C^{k+2})},$$

with $`G_k^{(i,j)}=G_k^{(j,i)}`$ for $`i>j`$.

We note two variants of (4.2), also valid for $`0ijk+1`$:

$`G_k^{(i,j)}`$ $`=x^{(ji)/2}(1+C)^{ji+1}{\displaystyle \frac{(1C^{i+1})(1C^{kj+1})}{(1C)(1C^{k+2})}},`$
$`=x^{1/2}C^{(ji+1)/2}{\displaystyle \frac{(1C^{i+1})(1C^{kj+1})}{(1C)(1C^{k+2})}}.`$
It is well known that $`G_k^{(i,j)}`$ is $`x^{(ji)/2}`$ times a rational function of $`y`$, and it is useful to have an explicit formula for it as a quotient of polynomials. (See Sato and Cong  and Krattenthaler .) Let us define polynomials $`p_n=p_n(x)`$ by

$$p_n(x)=\underset{0kn/2}{}(1)^k\left(\genfrac{}{}{0pt}{}{nk}{k}\right)x^k.$$

The first few values are

$`p_0`$ $`=1`$
$`p_1`$ $`=1`$
$`p_2`$ $`=1x`$
$`p_3`$ $`=12x`$
$`p_4`$ $`=13x+x^2`$
$`p_5`$ $`=14x+3x^2`$
$`p_6`$ $`=15x+6x^2x^3`$
These polynomials can be expressed in terms of the Chebyshev polynomials of the second kind $`U_n(x)`$ by

$$p_n(x)=x^{n/2}U_n\left(\frac{1}{2\sqrt{x}}\right).$$

It is not difficult to show that

$$p_n=\frac{1C^{n+1}}{(1C)(1+C)^n},$$

and thus we obtain

$$G_k^{(i,j)}=x^{(ji)/2}\frac{p_ip_{kj}}{p_{k+1}},$$

for $`0ijk`$, and in particular, $`G_k^{(j)}=x^{j/2}p_{kj}/p_{k+1}`$ and $`G_k=p_k/p_{k+1}.`$

We can now describe our combinatorial interpretation of $`T(3,n)`$: $`T(3,n)`$ counts triples of ballot paths whose heights are $`k`$, $`k2`$, and $`k4`$ for some $`k`$, ending at levels 4, 3, and 2, together with some additional paths of height at most 5. (Note that if a path of height $`k4`$ ends at level 2, then $`k`$ must be at least 6.) More precisely, let $`H_k^{(j)}`$ be the generating function for ballot paths of height $`k`$ that end at level $`j`$, so that $`H_k^{(j)}=G_k^{(j)}G_{k1}^{(j)}`$. Then we have:

###### Theorem 4.2.

(4.3)
$$1+\underset{n=0}{\overset{\mathrm{}}{}}T(3,n+1)x^n=\sqrt{x}\underset{k=6}{\overset{\mathrm{}}{}}H_k^{(4)}H_{k2}^{(3)}H_{k4}^{(2)}+2G_1+2G_2+G_3+G_5.$$

###### Proof.

First note that $`T(m,n)=\frac{1}{2}(1)^n4^{m+n}\left(\genfrac{}{}{0pt}{}{m{\scriptscriptstyle \frac{1}{2}}}{m+n}\right)`$, so

$`\frac{1}{2}(1)^m(14x)^{m1/2}`$ $`={\displaystyle \underset{n=m}{\overset{\mathrm{}}{}}}\frac{1}{2}(1)^n4^{m+n}\left({\displaystyle \genfrac{}{}{0pt}{}{m\frac{1}{2}}{m+n}}\right)x^{m+n}`$
$`={\displaystyle \underset{n=m}{\overset{0}{}}}\frac{1}{2}(1)^n4^{m+n}\left({\displaystyle \genfrac{}{}{0pt}{}{m\frac{1}{2}}{m+n}}\right)x^{m+n}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}T(m,n)x^{m+n}.`$
Setting $`m=3`$, dividing both sides by $`x^4`$, and subtracting the first sum on the right from both sides gives

(4.4)
$$\frac{(14x)^{5/2}}{2x^4}\frac{10}{x}+\frac{15}{x^2}\frac{5}{x^3}+\frac{1}{2x^4}=\underset{n=0}{\overset{\mathrm{}}{}}T(3,n+1)x^n.$$

Using the method described in Lemma 3.3 and Theorem 3.4, we find, with the help of Maple, that the sum

$$\sqrt{x}\underset{k=6}{\overset{\mathrm{}}{}}H_k^{(4)}H_{k2}^{(3)}H_{k4}^{(2)}$$

is equal to

$$\frac{(14x)^{5/2}}{2x^4}\frac{10}{x}+\frac{15}{x^2}\frac{5}{x^3}+\frac{1}{2x^4}+1\frac{2}{1x}2\frac{1x}{12x}\frac{12x}{13x+x^2}\frac{14x+3x^2}{15x+6x^2x^3}.$$

Then (4.3) follows from (4.4), the formula $`G_k=p_k/p_{k+1}`$, and the formulas for $`p_k`$, $`k=1,\mathrm{},6.`$


format: text

text:
# Directed-loop Monte Carlo simulations of vertex models

## I Introduction

Vertex models have a long and distinguished history in physics. Their fame is intimately connected to the concept of integrability, and the exact solutions of the six-vertex Lieb  and the eight-vertex Baxter-82  models with periodic boundary conditions (PBC) are indeed milestones in physics. Despite being exactly solvable, there are questions about these models that cannot easily be answered. An example is the influence of boundary conditions on correlation functions. While boundary conditions are not normally important in the thermodynamic limit, they have a profound influence on the vertex models. Exact studies, made for the six-vertex model with the domain wall boundary conditions (DWBC) KBI-93  show this in particular. These studies were restricted to certain points in the phase diagram, and involve rather sophisticated mathematical methods. It is thus appropriate to complement them with Monte Carlo simulations.

The purpose of this article is to demonstrate that the directed-loop Monte Carlo algorithm developed for quantum spin systems SS  can be used as an effective tool to study vertex models. The discussion of the algorithm will be kept general, but when demonstrating its use we will focus on the six-vertex model with the DWBC, a model which is difficult to simulate using other known Monte Carlo algorithms.

## II Monte Carlo algorithm

In a vertex model, each vertex have edges with an Ising-like variable, an arrow, that points either away from or into the vertex. The arrangement of arrows around the vertex determines the vertex weight. Two vertices are joined by their common edge, sharing the arrow on the edge. In general there are no restrictions on which vertices are joined, however for traditional vertex models nearest-neighbor vertices are joined together. The Monte Carlo algorithm discussed here always flips two (or zero) arrows on a vertex, thus it is limited to models where an even number of arrows are pointing away from each vertex. Most vertex models of interest obey this rule.

In visualizing the directed-loop Monte Carlo algorithm, originally developed for quantum systems SS , it is helpful to cut every edge into two pieces, each piece having an arrow belonging to a specific vertex, Fig. 1.

For a valid vertex configuration the arrows on the two parts of an edge must have the same orientation. The directed-loop algorithm is as follows: Pick a random vertex $`v_1`$ and a random edge belonging to that vertex. Based on these choices select in a probabilistic manner another edge belonging to $`v_1`$ and name that the out-edge. Then flip the arrows on both the part of the in-edge and the part of the out-edge belonging to $`v_1`$. This introduces two discontinuities in the arrow configurations on the edges, one on the starting in-edge and another one on the out-edge. The new configuration is thus not an allowed vertex configuration. To repair this, the out-edge discontinuity is moved by repeating the procedure on the vertex connected to the out-edge $`v_2`$, this time using the out-edge of $`v_1`$ and the in-edge on $`v_2`$. The process is stopped when the out-edge selected is the starting edge, thus healing all discontinuities. In this way arrows are flipped as a loop is constructed, and a new allowed vertex configuration is arrived at when the loop closes.

In order to determine the probabilities for selecting out-edges and to see how detailed balance is satisfied one needs to consider also the probability for the reverse update. The reverse update consists of traversing the same loop in the opposite direction while flipping arrows back. As is explained in detail in Ref. SS , detailed balance is satisfied for the whole loop construction, if detailed balance is satisfied in each edge selecting step, for which the criterion is as follows: Let $`w`$ be the weight of the vertex $`v`$ before edge-flips, then the probability $`P(v,io)`$ for exiting at the out-edge $`o`$, given that the in-edge is $`i`$, should satisfy

$$wP(v,io)=w^{}P(v^{},oi),$$
(1)
where $`w^{}`$ is the weight of the vertex $`v^{}`$ obtained by flipping the arrows on edges $`i`$ and $`o`$ belonging to the vertex $`v`$. Notice that $`P(v^{},oi)`$, on the right hand side, describes an edge-selecting step in the reverse update process where the loop is traversed in the opposite direction to that described on the left hand side. The loop construction should not terminate in the edge-selecting step, thus

$$\underset{o}{}P(v,io)=1,$$
(2)
where the sum is taken over all possible out-edges, including the in-edge $`i`$.

This algorithm resembles closely the ice model algorithm invented by Rahman and Stillinger Rahman , generalized to arbitrary couplings by Barkema and Newman BN . In fact, at the point in parameter space where all vertex weights are equal our algorithm is identical to the long-loop version of the ice model algorithm. However away from this point, Barkema and Newman’s algorithm involves accepting or rejecting the loop after it has been constructed. The directed-loop algorithm has no such accept/reject step. A comparison of integrated autocorrelation times for the directed-loop algorithm and the short-loop algorithm of Barkema and Newman are shown in Fig. 3. The autocorrelation times are measured in units of lattice sweeps. One lattice sweep corresponds to a number of elementary loop moves such that on average each vertex on the lattice have been visited twice. In defining visited we include parts of the loop where the loop bounces off a vertex (relevant for the directed-loop algorithm) and the neck part of short-loops. Neither the bounces nor the short-loop-necks contribute to changes in the vertex configuration. However they are intrinsic parts of the algorithms and requires computer effort, and should therefore be accounted for.

The upper panel of Fig. 3 shows integrated autocorrelation times of the observable counting the number of $`c`$-type vertices in each configuration. This observables was chosen to compare with the performance results in Ref. BN . While the integrated autocorrelation times are larger for the short-loop algorithm the scaling with system size appears to be equal for both algorithms. The lower panel shows integrated autocorrelation times for the total arrow-polarization in the y-direction. These scales much worse for the short-loop algorithm than for the directed-loop algorithm. This is to be expected from the fact that most loops accepted in the short-loop algorithm are small, while large loops that wind around the boundary of the lattice is needed to change the total polarization. These are not suppressed in the directed-loop algorithm, thus leading to better performance.

The Eqs. (1) and (2) form several coupled sets which in most cases are under-determined. There are thus many solutions for the out-edge selection probabilities $`P`$. Some general solutions and analysis of their efficiency for different quantum systems were reported in Ref. Olav . Here we employ the solution B in Ref. Olav  to the eight-vertex model, but solutions for higher-vertex models are not hard to find as well. The allowed vertices for the eight-vertex model and their statistical weights are shown in Fig. 2. To shorten notation, we consider the so-called symmetric case: the statistical weights, $`a,`$ $`b,`$ $`c,`$ and $`d,`$ of the allowed states are assumed to be invariant under the simultaneous reversal of all arrows. The generalization of the algorithm to the non-symmetric case can be performed easily.

Let $`W_1,\mathrm{},W_4`$ be the vertex weights $`a,b,c,d`$ of the eight-vertex model ordered so that $`W_1W_2W_3W_4`$. Then the probability for picking the out-edge on a vertex with weight $`W_i`$ resulting in a new vertex weight $`W_j`$ after flipping arrows is $`t_{ij}/W_i,`$ where $`t_{ij}=t_{ji}`$ and the non-zero entries of the $`4\times 4`$ matrix $`t`$ are

$`t_{12}`$ $`=`$ $`(W_1+W_2W_3W_4)/2,`$
$`t_{13}`$ $`=`$ $`(W_1W_2+W_3W_4)/2,`$ (3)
$`t_{23}`$ $`=`$ $`(W_1+W_2+W_3+W_4)/2,`$
$`t_{14}`$ $`=`$ $`W_4,`$
when $`W_1W_2W_3W_40`$. Otherwise one needs to include bounces in which the out-edge coincides with the in-edge. In this case a solution can be chosen as follows:

$`t_{11}`$ $`=`$ $`W_1W_2W_3W_4,`$
$`t_{1j}`$ $`=`$ $`t_{j1}=W_j,j=2,3,4,`$ (4)
$`t_{ij}`$ $`=`$ $`0,\mathrm{otherwise}.`$
The directed-loop algorithm satisfies ergodicity as any configuration can be obtained from another configuration by flipping spins along a finite number of (possibly overlapping) loops. This process is exactly the directed-loop update, thus ergodicity follows.

The algorithm presented here has many similar features to the loop algorithm Evertz . The loop algorithm breakup rules for the symmetric eight-vertex model can be chosen identical to Eq. (3), as can be seen from Ref. Kawashima , whenever the weights are such that no bounces are needed in the directed-loop algorithm. However in parameter regimes where bounces are needed, the related feature in the loop algorithm is to “freeze” independent loops together. Bounces and “freezing” of loops are very different in how they act to change the configuration. While bounces is a local resistance to changing a vertex, “freezing” causes big non-local changes of the vertex configuration. There are also other differences: For general vertex models the set of non-freezing/bouncefree solutions is always smaller for the loop algorithm than for the directed-loop algorithm. This allows the directed-loop algorithm to be efficient in a larger region of parameter space than the loop algorithm. In particular this applies to the asymmetric eight-vertex model.

Note that the need for bounces is generally not so crucial for higher-vertex models with many weights of the same magnitude, thus we expect that the directed-loop algorithm should work well in simulating these. Note also that an algorithm based on the directed-loop idea was recently demonstrated to be effective in simulating classical integer-valued link-current models Alet .

## III Six-vertex model with the DWBC

The six-vertex model with the DWBC was introduced in Ref. K-82  in connection with the calculation of the correlation functions for exactly solvable $`1+1`$ dimensional models KBI-93 .

Here we recall the definition of the model in brief, referring for further details to the Ref. BPZ-02  where a more detailed description of the model and a comprehensive list of references are given.

The model is defined on an $`N\times N`$ square lattice; the thermodynamic limit corresponds to $`N\mathrm{}.`$ There are six possible states at each vertex: one should set $`d=0`$ in the eight-vertex model defined above, Fig. 2. The model is symmetric: the statistical weights, $`a,`$ $`b,`$ and $`c,`$ of the allowed states are assumed to be invariant under the simultaneous reversal of all arrows. Hence, the model is characterized by only two parameters, which can be taken to be $`a/c`$ and $`b/c`$. We set $`c=1`$ henceforth.

The DWBC imply that all arrows on the top and bottom of the lattice are pointing inward, while all arrows on the left and right boundaries are pointing outward, Fig. 4.

To investigate the spatially inhomogeneous behavior of this model we focus on the polarization, $`\chi _N(x,y)`$ BPZ-02 ; BKZ-02 , which is the ensemble average of the arrow direction on the edge with coordinates $`(x,y)`$ on the $`N\times N`$ lattice. The coordinate system used is shown in Fig. 4. Due to the symmetry of the model it is sufficient to consider the polarization of the horizontal arrows only. The value $`+1`$ ($`1`$) is assigned to an arrow pointing to the right (left) and the ensemble average is assumed to be normalized by dividing by the partition function. Therefore, $`\chi _N`$ lies between $`1`$ and $`1.`$

Obviously, $`\chi _N`$ is independent of the coordinates of the edge in case of PBC. For these boundary conditions $`\chi _N`$ is known in the thermodynamic limit, and exhibits ferroelectric order, antiferroelectric (AF) order or no order, depending on the position on the $`(a,b)`$ plane. Thus, three phases exist in the six-vertex model with PBC: ferroelectric, antiferroelectric, and disordered phase. In Fig. 5 the phase diagram on the $`(a,b)`$ plane for the model with PBC is plotted (cf., Fig. 8.5 of Ref. Baxter-82 ).

Introduce a parameter $`\mathrm{\Delta }`$ by the formula

$$\mathrm{\Delta }=\frac{a^2+b^21}{2ab}.$$
(5)
The case $`\mathrm{\Delta }>1`$ (regions $`\mathrm{I}`$ and $`\mathrm{II}`$ in Fig. 5) corresponds to the ferroelectric phase; the case $`1<\mathrm{\Delta }<1`$ (region $`\mathrm{III}`$ in Fig. 5) to the disordered phase; the case $`\mathrm{\Delta }<1`$ (region $`\mathrm{IV}`$ in Fig. 5) to the AF phase.

Fig. 5 may be considered as the phase diagram for the model with the DWBC, in the sense that the free energy takes a different analytic form in the regions $`\mathrm{I}`$ through $`\mathrm{IV}`$ (see Ref. ZJ-00  for details). But, in case of the DWBC the polarization $`\chi _N`$ depends on the position of the edge. In the next section we show numerical results for the polarization $`\chi _N(x,y)`$ of the horizontal arrows as the parameters $`a`$ and $`b`$ are varied.

Making use of the directed-loop algorithm described in the previous section for simulation of the model with the DWBC one should treat vertices belonging to the boundary and the bulk vertices differently. In the bulk one finds by setting $`d=W_4=0`$ and $`c=1`$ in Eqs. (3), that bounces are only necessary when $`a+b<1`$ or $`|ab|>1`$. For the boundary vertices the loop is not allowed to exit on the boundary edges, because the arrows on these edges are fixed by the boundary conditions. This leads to more restricted equation sets (many $`W`$’s are equal to zero) for the boundary vertices and generally requires the inclusion of bounce processes.

Another important point should be mentioned is that the DWBC do not violate the ergodicity of the algorithm even though loops which wind around the boundaries are excluded. These winding loops are needed in order to change the net polarization in the $`x`$\- or $`y`$-direction. However, one can verify that the boundary conditions restricts the net polarization in both these direction to be zero for any configuration, so winding loops are not necessary to sample the full configuration space allowed by the boundary conditions.

## IV Results

In this section we present the results of the simulations for the polarization $`\chi _N(x,y)`$ in the disordered, antiferroelectric and ferroelectric phases.

(i) Disordered phase: $`1<\mathrm{\Delta }<1.`$ First consider the particular case $`\mathrm{\Delta }=0`$ (dotted quartercircle in Fig. 5). An exact expression for $`\chi _N(x,y)`$ in this case was obtained by Kapitonov and Pronko KP-04  recently. To check our algorithm we have compared results for the polarization at the point $`a=b=1/\sqrt{2}`$ with the exact results of Ref. KP-04 . The comparison can be seen in Fig. 6,

where the polarization is shown as a function of $`x`$ for different values of $`y`$ and system sizes, $`N`$. One can clearly see that the boundary values of the polarization, $`\pm 1,`$ extends a finite distance into the bulk and this distance depends on $`y.`$ The areas where the polarization stays at its boundary values are termed “frozen” regions. Going further into the bulk, there is a transition to a “disordered” region, where apart from small wiggles due to the finite system size, the polarization changes smoothly. It is interesting to note that there never is any extended regime where the polarization is zero, as is the case for PBC. The transition between the “frozen” and “disordered” regions gets sharper as the system size is increased, as can be seen by comparing the two panels in Fig. 6.

It is convenient to visualize the behavior of the polarization using greyscale plots, where greyvalues are assigned to values of $`\chi _N(x,y)`$ and each point $`(x,y)`$ corresponds to a location of the midpoint of a horizontal edge following the layout described in Fig. 4. For $`a=b=1/\sqrt{2}`$ such a plot is shown in Fig. 7(a).

The four “frozen” corners are clearly apparent. In these regions, the vertices are all of the same type, and are, from upper left to bottom right, $`a_1`$, $`b_1`$, $`b_2`$, $`a_2`$, respectively. To measure the area of the “frozen” regions, we define a threshold value $`ϵ=0.08`$, such that points $`(x,y)`$ where $`|\chi _N(x,y)|>1ϵ`$ are judged to be in a “frozen” region. Applying this we find that each “frozen” corner is 4.6% of the total area. This value changes relatively little changing the value of $`ϵ`$.

Going away from the $`\mathrm{\Delta }=0`$ curve, let us follow along the diagonal, $`a=b,`$ towards $`\mathrm{\Delta }=\mathrm{}`$ first, Fig. 7. As the values of the vertex weights $`a`$ and $`b`$ increase, the area of the “frozen” regions decreases. We find that with $`ϵ=0.08`$ each frozen corner in $`(b)`$ is $`4.0\%`$ of the total area, and $`2.8\%`$ in (c). For very large values of $`a=b`$, the polarization $`\chi _N(x,y)`$ increases linearly from $`1`$ to $`1`$ as $`(x1)/N`$ goes from $`0`$ to $`1`$, independent of $`y`$, as can be seen in Fig. 7 (d). This is consistent with what is expected from an ensemble of configurations with the smallest possible number of $`c`$-type vertices: $`N!`$ configurations each with a single $`c`$-type vertex on every row and column.

Consider now $`ab.`$ Because of the symmetry of the phase diagram, Fig. 5, one can choose $`b>a`$ without loss of generality. The weights of the vertices in the four “frozen” corners are no longer equal, and the “disordered” region distorts into an oblong shape oriented along the diagonal with large corners of $`b_2`$ and $`b_1`$ vertices, see Fig. 8. The simulations for $`a=1/4`$ and $`b=\sqrt{15/16}`$ are shown in Fig. 8(a). The width of the oblong region shrinks as $`b`$ increases keeping $`a`$ fixed, $`a=1/4`$, and becomes very thin at the boundary to the ferroelectric region, as can be seen in Fig. 8(b). Along this boundary, $`b=a+1`$, the width of the oblong region expands as $`a`$ increases with $`N`$ being constant.

(ii) Antiferroelectric phase: $`\mathrm{\Delta }<1.`$ The simulations in the AF phase are less efficient than in the disordered phase. This is partly due to the presence of the bounce processes also for bulk vertices, but another feature which makes the simulations difficult in this phase is the degeneracy of the two types of AF orders. In the AF phase it becomes energetically favorable to have a maximum possible amount of $`c`$-type vertices, which is achieved by placing $`c`$-type vertices in a diamond placed in the center of the lattice. For an even $`N`$ this diamond can be placed in two equivalent places differing only by one lattice spacing, as shown in Fig. 9.

The Monte Carlo algorithm is however slow in tunnelling between these configurations, and this sets a limit to its performance. For odd $`N`$ there is no such a degeneracy and the simulations are more efficient. Greyscale plots of the polarization for $`a=b=1/2`$ and $`a=b=3/8`$ are shown in Fig. 10.

We have plotted results for both even and odd $`N.`$

One can see that the “disordered” region have a diamond-like shape, which is consistent with the domination of the $`c`$-type vertices in the AF phase. As $`a=b`$ decreases $`(\mathrm{\Delta }\mathrm{})`$, the shape of the “disordered” region should converge to the one shown in Fig. 9, that is, the boundaries of the “disordered” region should become more and more straight. But, this convergence appears to be rather slow and it is not easy to see it from Fig. 10. What one can clearly see from Fig. 10 is the difference between odd and even $`N.`$ For odd $`N`$ AF oscillations are clearly visible in the center of Figs. 10(c) and (d), while they are much weaker for even $`N`$, Figs. 10(a) and (b), reflecting the degeneracy mentioned above. These difference between even and odd $`N`$ can also be clearly seen from Fig. 14. For odd $`N`$ AF oscillations are weaker at $`a=b=1/2`$ than at $`a=b=3/8`$.

For $`ab`$ greyscale plots are shown in Fig. 11.

Here AF oscillations in the middle of the plot are visible for $`a=1/4`$ and $`b=1/2`$, Fig. 11(a), while they have almost vanished at the boundary between the AF and disordered phases, Fig. 11(b).

(iii) Ferroelectric phase: $`\mathrm{\Delta }>1.`$ The behavior of the polarization in this phase is essentially the same as shown in Fig. 8. Vertices of type $`b`$ dominate completely in the region $`\mathrm{II}`$ of the phase plane Fig. 5, while in the region $`\mathrm{I}`$ of the phase plane the dominant vertices are those of type $`a.`$ If one goes along the phase boundary, $`b=a+1,`$ towards $`a=\mathrm{},`$ the widths of the “disordered” region is increased, as we have mentioned in the end of the part (i) of this Section.

The exact expression is known BPZ-02  for the polarization along the boundary, $`\chi _N(x,1).`$ Comparing our Monte Carlo data to this expression we find that in no cases is the absolute difference bigger than $`0.016,`$ which is comparable to the statistical errors of our simulations.

## V Discussion

We have considered the phase diagram of the model for the given $`N.`$ Now, discuss the following problem: what happens with $`\chi _N(x,y)`$ in the thermodynamic limit, $`N\mathrm{}\mathrm{?}`$ It is natural to expect that differences in the behavior of the polarization in the different parts of the phase plane, Fig 5, become more pronounced as $`N\mathrm{}.`$ As one can see in Fig. 6, the wiggles in the “disordered” region decrease with $`N`$ increasing, and this is, indeed, the case for all the points ($`a,b`$) lying in the disordered phase ($`1<\mathrm{\Delta }<1,`$ region $`\mathrm{III}`$ of the phase plane, Fig. 5) and checked in our simulations. We expect that these wiggles, coming from the antiferroelectrically ordered configurations, should vanish completely in the thermodynamic limit for this phase. The next conjecture we want to make is on the behavior of the polarization along the boundary, $`\chi _N(x,1).`$ It is known that for $`\mathrm{\Delta }=0,`$ as well as at the point $`a=b=1,`$ the boundary polarization becomes the Heaviside step function in the thermodynamic limit BPZ-02 ; Z-96 . We conjecture that this is the case for the whole disordered phase; the position of the discontinuity will depend on the ratio between $`a`$ and $`b.`$ We present Fig. 12 to support this conjecture.

Furthermore, note that for $`a=b=1/\sqrt{2}`$ there is a mapping (see, e.g., Ref. ZJ-00 ) of the six-vertex model with the DWBC onto the so-called model of domino tilings of the Aztec diamond. The thermodynamic behavior of the latter model was investigated in Refs. JPS-98 . It shows the same features as in Fig. 7(a): the tilings are ordered (frozen) in the corners of the diamond, while going inside one falls into the “disordered” region. All these features were named the “Arctic Circle Theorem”, since the shape of the boundary between the “frozen” and “disordered” regions is circular. The transition between “frozen” and “disordered” regions is step-like, with the height of the step function depending on the coordinates $`x`$ and $`y.`$

We expect the analogue of the Arctic Circle Theorem to take place for the whole disordered phase, $`1<\mathrm{\Delta }<1`$: there should be the “frozen” regions, “disordered” region, and a sharp transition between them. We expect also that the profile of the boundary between the “frozen” and “disordered” regions is circular for $`a=b`$, even though there is no obvious symmetry protecting this statement. Note that the very “smeared” profile in Fig. 7(d) does not contradict our hypotheses because $`N=64`$ is relatively small compared to the values of the vertex weights $`a`$ and $`b,`$ and is thus far from the thermodynamic limit for this point of the phase diagram.

For the ferroelectric phase, $`\mathrm{\Delta }>1,`$ the greyscale plot Fig. 8(b) together with the scans shown in Fig. 13 leads to the natural conjecture: in the whole region $`\mathrm{II}`$ of the phase plane, Fig. 5, a sharp discontinuity from a “frozen” domain with $`b_1`$ vertices to the one with $`b_2`$-vertices takes place in the thermodynamical limit. In the region $`\mathrm{I}`$ the behavior is essentially the same, one should simply use $`a`$-type vertices instead of the $`b`$-type.

To this end, consider the antiferroelectric phase, $`\mathrm{\Delta }<1`$. We expect the step-like behavior of the boundary polarization, $`\chi _N(x,1)`$ in this phase in the thermodynamic limit, as well as the existence of the “frozen” regions in the corners. Our statements on the behavior of the polarization deep inside the lattice are more speculative. For $`a=b`$ and even $`N`$ the height of the AF oscillations decreases, while for odd $`N`$ these oscillations seem to be non-vanishing in the large $`N`$ limit, see Fig. 14.

Our belief is that there is a finite region with AF order for odd $`N,`$ as $`N\mathrm{},`$ while for even $`N`$ the polarization exhibits no such an order.

Finally, we would like to stress that the directed-loop algorithm can also be applied to study the six-vertex model with any boundary conditions, and the higher-vertex models. These could help in solving the problems for which the analytical methods are difficult to apply. For example, the six-vertex model with any boundary conditions can be considered as a model for a description of interface roughening of a crystal surface Beij-77 . An important point in studies Refs. Beij-77  is the existence of exact analytical results for the six-vertex model with PBC Lieb ; Baxter-82 . Therefore, numerical data referring to other boundary conditions than PBC could give a new insight for these studies.

###### Acknowledgements.

We thank V.V. Cheianov and A.G. Pronko for useful discussions and the authors of the work KP-04  for providing us with their data. M.B. Zvonarev’s work was supported by the Danish Technical Research Council via the Framework Programme on Superconductivity. Monte Carlo calculations were in part carried out using NorduGrid, a Nordic facility for Wide Area Computing and Data Handling.


text:
# A unique representation of polyhedral types. Centering via Möbius transformations

## Proof of Lemma 1.

The Möbius transformations of $`S^d`$ correspond to isometries of the hyperbolic space $`H^{d+1}`$, of which $`S^d`$ is the infinite boundary.

For $`n3`$ points $`v_1,\mathrm{},v_nS^d`$, we are going show that there is a unique point $`xH^{d+1}`$ such that the sum of the ‘distances’ to $`v_1,\mathrm{},v_n`$ is minimal. Of course, the distance to a point in the infinite boundary is infinite. The right quantity to use instead is the distance to a horosphere through the infinite point (see the figure).

###### Definition.

For a horosphere $`h`$ in $`H^{d+1}`$, define

$$\delta _h:H^{d+1},$$

$$\delta _h(x)=\{\begin{array}{cc}\mathrm{dist}(x,h)\hfill & \text{if }x\text{ is inside }h,\hfill \\ 0\hfill & \text{if }xh,\hfill \\ \mathrm{dist}(x,h)\hfill & \text{if }x\text{ is outside }h,\hfill \end{array}$$

where $`\mathrm{dist}(x,h)`$ is the distance from the point $`x`$ to the horosphere $`h`$.

Suppose $`v`$ is the infinite point of the horosphere $`h`$. Then the shortest path from $`x`$ to $`h`$ lies on the geodesic connecting $`x`$ and $`v`$. If $`h^{}`$ is another horosphere through $`v`$, then $`\delta _h\delta _h^{}`$ is constant. If $`g:H^{d+1}`$ is an arc-length parametrized geodesic, then $`\delta _hg`$ is a strictly convex function, unless $`v`$ is an infinite endpoint of the geodesic $`g`$. In that case, $`\delta _hg(s)=\pm (ss_0)`$. These claims are straightforward to prove using the Poincaré half-space model, where hyperbolic space is identified with the upper half space:

$$H^{d+1}=\left\{(x_0,\mathrm{},x_d)^{d+1}\right|x_0>0\},$$

and the metric is

$$ds^2=\frac{1}{x_0^2}\left(dx_0^2+dx_1^2+\mathrm{}+dx_d^2\right)$$

Also, one finds that, as $`xH^{d+1}`$ approaches the infinite boundary,

$$\underset{x\mathrm{}}{lim}\underset{j=1}{\overset{n}{}}\delta _{h_j}(x)=\mathrm{},$$

where $`h_j`$ are horospheres through different infinite points and $`n3`$. Thus, the following definition of the point of minimal distance sum is proper.

###### Lemma (and Definition) 2.

Let $`v_1,\mathrm{},v_n`$ be $`n`$ points in the infinite boundary of $`H^{d+1}`$, where $`n3`$. Choose horospheres $`h_1,\mathrm{},h_n`$ through $`v_1,\mathrm{},v_n`$, respectively. Then there is a unique point $`xH^{d+1}`$ for which $`_{j=1}^n\delta _{h_j}(x)`$ is minimal. This point $`x`$ does not depend on the choice of horospheres. It is the point of minimal distance sum from the infinite points $`v_1,\mathrm{},v_n`$.

In the Poincaré ball model, hyperbolic space is identified with the unit ball as in the Klein model, but the metric is $`ds^2=\frac{4}{(1{\scriptscriptstyle x_j^2})_{\text{}}^2}dx_j^2.`$ (Since the Klein model and the Poincaré ball model agree on the infinite boundary and in the center of the sphere, one might as well use the Klein model in the following lemma.)

###### Lemma 3.

Let $`v_1,\mathrm{},v_n`$ be $`n3`$ different points in the infinite boundary of $`H^{d+1}`$. In the Poincaré ball model, $`v_jS^d^{d+1}`$. The origin is the point of minimal distance sum, if and only if $`v_j=0`$.

###### Proof.

If $`h_j`$ is a horosphere through $`v_j`$, then the gradient of $`\delta _{h_j}`$ at the origin is the unit vector $`\frac{1}{2}v_j`$. ∎

Lemma 1 is now almost immediate. Let $`x`$ be the point of minimal distance sum from the $`v_1,\mathrm{},v_n`$ in the Poincaré ball model. There is a hyperbolic isometry $`T`$ which moves $`x`$ into the origin. If $`\stackrel{~}{T}`$ is another hyperbolic isometry which moves $`x`$ into the origin, then $`\stackrel{~}{T}=RT`$, with $`R`$ is an orthogonal transformation of $`^{d+1}`$.

This concludes the proof of Lemma 1.

## Acknowledgements

I would like to thank Alexander Bobenko and Günter Ziegler for making me familiar with the problem of finding unique representatives for polyhedral types, and Ulrich Pinkall, who has provided the essential insight for this solution.


text:
# Statistical Challenges of Cosmic Microwave Background Analysis

## I Overview

What is cosmic microwave background (CMB) statistics and what is challenging about it?222For online material relating to this talk please refer to
www-conf.slac.stanford.edu/phystat2003/talks/wandelt/invited/ It involves estimating the covariance structure of of a spatial random field with $`10^6`$$`10^8`$ pixels, given only ONE realization of this field. The covariance matrix of these pixels is a complicated non-linear function of the physical parameters of interest. Of these physical parameters there are between 10 and 20, so even finding the maximum likelihood point is hard—determining and summarizing confidence intervals around the maximum likelihood point is very non-trivial. Cosmologists want to do all this and have the option of building in exact or approximate physical constraints on relationships between parameters. In addition, since collecting cosmological data is so difficult and expensive we want to combine all available data sets—both to test them for mutual disagreement which might signal new physics, and to improve the parameter inferences. In all of this the quantification of the uncertainties in the results is extremely important—after all the stated significance of our results will either drive or stop theoretical investigations and the design of new observational campaigns.

Before I get on to CMB specifics in section II let me give you the short version of (most) of this talk for statisticians: “The CMB is an isotropic (homoschedastic) Gaussian random field $`s`$ on the sphere. The desired set of cosmological parameters $`\mathrm{\Theta }=\{\theta _ii=1,\mathrm{},n\}`$, are related in a non-linear way to the spatial covariance structure $`S_{ij}s_is_j`$ of the field. Observers present us with a sampled, noisy, filtered and censored/polluted measurement of this field in several ’colors’. The analysis task is two-fold: infer the covariance structure of the field $`s`$. Infer the parameters $`\mathrm{\Theta }`$.” This is what could be termed “canonical” CMB analysis.

In this talk I will mainly describe challenges presented by this canonical CMB analysis. After a brief review of the scientific motivation for studying the CMB in section II I will describe the form of CMB data as well as the current status and prospects of obtaining it in section III. Section IV then outlines a framework for extracting cosmologically useful information from the data and section V illuminates some examples of challenges that arise when implementing this framework. I will touch on statistical questions concerning “non-canonical” CMB analysis in section VI and then conclude in section VII.

So why are we interested in facing the statistical challenges of CMB analysis?

## II What can we learn from the CMB?

Cosmologists are interested in studying the origins of the physical Universe. In order to do so they have to rely on data. For cosmologists, one of the great practical advantages of Einstein’s relativity over Newtonian physics is the fact that we cannot help but look into the past. Therefore, by observing light that reaches us from farther and farther away, we can study the Universe directly at earlier and earlier times, at least to the extent to which the Universe is transparent to light. Since the early Universe was a hot and opaque plasma we can only see back to the time when the plasma cooled sufficiently (due to the Hubble expansion) to combine into neutral atoms and the mean free time between photons collisions became of order of the present age of the Universe. Photons that we observe today which scattered for the last time in the primordial plasma are the CMB.

The change from plasma to gas, happened when the Universe was approximately 380,000 years old. CMB photons emitted at this time are therefore the most direct messengers that we can detect today of the conditions present in the Universe shortly after the Big Bang.333There are other messengers, namely neutrinos and gravitational waves, reaching us from even earlier times but we do not (yet) have the technology to detect them in relevant quantities. They constitute a pristine snapshot of the infant Universe which provide us with direct cosmological information uncluttered by the complex non-linear physics which led to the formation of stars and galaxies. One of the main attractions of the CMB is the conceptual simplicity with which it can be linked to the global properties of the Universe and the physics which shaped it at or near the Planck scale.

As a simple example, the serendipitous discovery of the CMB by Penzias and Wilson in 1965 showed that the CMB is isotropic to a very high degree. This was one of the key motivations for the development of the inflationary paradigm Guth ; Linde ; AlbrechtSteinhardt . Inflation describes generically the emergence (from the era of quantum gravity) of a large homogeneous and isotropic Universe. Inflation also predicted the spectrum of small metric perturbations from which later structure developed through the gravitational instability (though it was not the only mechanism to do so). The corresponding anisotropies in the CMB were first convincingly detected by the DMR instrument on the COBE satellite in 1992. The rejection of alternative mechanisms for the generation of the primordial spectrum of metric perturbations in favor of inflation was a major advance driven by the measurement of the large angle CMB anisotropy. These developments have led to inflation becoming part of the current cosmological standard model. A very robust prediction of generic model implementations of inflation is the Gaussianity and homogeneity of the resulting perturbations.

Quantitatively, the properties of our Universe are encoded in a set of $`n1020`$ cosmological parameters $`\mathrm{\Theta }`$ where $`n`$ depends on the level of detail of the modeling or, commonly, on the specification of theoretical priors which fix some of these parameters to “reasonable” values. These parameters specify the geometry and average energy density of the Universe, as well as the relative amounts of energy density contributed by the ingredients of the primordial soup (dark matter, ordinary baryons, neutrinos and dark energy and photons). In addition, the anisotropy carries information about the spectrum of primordial (inflationary) perturbations as well as their type (adiabatic or isocurvature). By combining observations of the anisotropy of both the effective temperature and the polarization of the CMB photons we can infer how transparent the Universe really was for the CMB photons on their way from last scattering to hitting our detectors. This in turn can tell us about the history of star formation.

A very exciting prospect is that by studying the details of CMB polarization we can infer the presence or absence of gravitational waves at the time of last scattering. A detection would offer an indirect view of one of the elusive messengers that started their journey at an even earlier epoch, adding a nearly independent constraint on the properties of the Universe at the Planck scale.

All this information is not encoded in actual features in the CMB map of (temperature or polarization) anisotropies. In fact in a globally isotropic universe the absolute placement of individual hot and cold spots is devoid of useful information. Information can, however, be stored in the invariants of the photon brightness fluctuations under the group of rotations SO(3). These are the properties of the field that only depend on the relative angular distance between two points of the field. For a Gaussian field, where 2-point statistics specify all higher order moments, this means that the angular power spectrum coefficients of the anisotropies contain all of the information.

The challenge for theoreticians was then to develop a detailed theory of the angular power spectrum $`C_{\mathrm{}}`$, as a function of angular wavenumber $`\mathrm{}`$, given the cosmological parameters $`\mathrm{\Theta }`$. While conceptually simple, it required a decade-long intellectual effort to model the relevant physical processes at the required level of precision. As a result, there now exist several Boltzmann codes (e.g. CMBFAST CMBFAST  or CAMB CAMB ), which numerically compute $`C_{\mathrm{}}(\mathrm{\Theta })`$ to 1% precision or better. The power spectra $`C_{\mathrm{}}(\mathrm{\Theta })`$ are sensitive functions of certain combinations of the parameters and weak functions of others (degeneracies). These weakly constrained parameter combinations are referred to as degeneracies. Within the context of the standard cosmological model, the theory of this dependence is well-understood.

It is clear from the preceding discussion that the CMB is an extremely valuable source of what amounts to “cosmological gold”: information about the physics and the global properties of the early Universe. So what are the observational prospects?

## III Data from CMB Observations

A major international effort is underway to make high quality observations of the microwave sky using ground-based, balloon borne and space missions CMBobs . Space missions have the advantage of being able to scan the whole microwave sky. NASA’s “Wilkinson Microwave Anisotropy Probe” (WMAP) was launched in 2001 with great success and is currently in operation. It reported its first year of data earlier this year WMAP . WMAP will continue to collect data for at least another three years. In the medium term (ie. in late 2007) we anticipate the launch of “Planck,” a joint ESA/NASA space mission.444http://www.rssd.esa.int/index.php?project=PLANCK Additional space missions may follow that will focus on measuring the polarization anisotropies in the CMB. In the meantime ground and balloon-based missions are jostling to accelerate our learning curve by providing maps at high angular resolution on small patches of the sky (up to a few degrees large). As a result, by the end of this decade we will have a mountain of CMB data.

From this mountain (and of course from the part of it which we have already available today) we would like to extract the cosmological gold. In order to do so we need to understand how the data and the information are related. In Figure 1 I show a simplified schematic of this relationship and I will now go through the various steps.

Starting with a set of cosmological parameters $`\mathrm{\Theta }`$ we can use a Boltzmann code to compute the power spectrum $`C_{\mathrm{}}`$. Given this power spectrum and the assumption of Gaussianity and isotropy we have all the information we need to create a statistical realization of a CMB anisotropy map. This is most simply done working in the Fourier representation: we draw the spherical harmonic coefficients $`a_\mathrm{}m`$ from uncorrelated Gaussians with zero mean and variance given by $`C_{\mathrm{}}`$. Then we compute the spherical harmonic transform to obtain

$$T(\widehat{n})=\underset{\mathrm{}=2}{\overset{\mathrm{}=\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{m=\mathrm{}}{}}T_\mathrm{}mY_\mathrm{}m(\widehat{n}),$$
(1)
where $`Y_\mathrm{}m(\widehat{n})`$ are the spherical harmonics, an orthonormal and complete basis for functions on the sphere. The unit vector $`\widehat{n}`$ points in a direction on the sky.

Unfortunately, the $`T_\mathrm{}m`$ or $`T(\widehat{n})`$ cannot be observed directly. We are tied to our location in the Galaxy. There are various sources of copious amounts of microwave radiation in the Galaxy, such as dust and electrons from winds which are accelerated in the Galactic magnetic field. These foregrounds add to the CMB signal to make our sky.

This sky is then observed in various frequency bands (indicated by the red, green and blue arrows) by a CMB instrument (the figure shows an artist’s conception of the WMAP satellite wmapsite ). This instrument itself has a complicated transfer function: since weight and size constraints force satellite optics to be built compactly the optics are not free from distortion. Microwaves have macroscopic wavelengths and therefore diffract around the edges of the instrument. This leads to sidelobes in the beam maps. The instrument scans the sky in a certain pattern (the “scan strategy”) and internal instrument systematics are added to the scanned signal. The microwave detectors (either radiometers or bolometers) add noise to generate the time ordered data (TOD).

There are different levels of detail of the CMB data. To give an idea for the orders of magnitude involved, let us go through the sizes of the various objects in Figure 1 for the Planck experiment.

* The complete TOD for Planck will take up of order 1 Terabyte (=$`10^{12}`$ bytes) of storage (without counting house-keeping and pointing data).
* Each of the 100 detectors (channels) results in a map which is of order 10-100 Megabytes.
* The channels are grouped into of order 10 frequency bands.
* The combined maps at these different frequency bands will be combined into maps of the physical components (such as dust, synchrotron, CMB).
* The CMB power spectrum has a few thousand coefficients $`C_{\mathrm{}}`$.
* These power spectrum coefficients are a function of $`1020`$ cosmological parameters.

Note that there is a trade-off between the level of compression and what assumptions are implemented in the data analysis. Note also, that except for the raw data each of these data products by themselves mean little unless some means of assessing their statistical uncertainty is provided.

We immediately run into practical problems. For example, if we would like to specify a noise covariance matrix for each combined map at each frequency we would have to specify 10 times $`(10^6)^2/2`$ elements of a matrix. The necessary storage space of order 10,000 Gigabyte basically precludes practical public distribution of the data. Other ways of specifying the uncertainty must be found.

Before we discuss in detail how cosmological information is encoded in the data, let us comment about an aspect of CMB statistics that is challenging even before there is any data in hand. This aspect is experimental design. While our experimental colleagues are putting a great deal of valuable thought into their instruments and observational strategies, it is currently still done in an informal way, by ingenuity rather than by formal method. This is partially so because it is hard to define optimality criteria that everyone would agree with and partially because of the immense effort involved to actually carry out the optimization. Evaluating any one proposed design requires many simulations of the full process from observation to analysis—the very process that presents the challenges I am discussing in this talk. Doing this repeatedly to search the space of design parameters for an optimal solution would be an immense task. Asymptotic techniques have been implemented on the basis of the Fisher matrix minimum variance predictions, but it is important to keep in mind that these are lower bounds on the expected variance, assuming a unimodal, Gaussian likelihood shape.

## IV Inference From the Data

A stochastic model of the data and the information contained in it can be summarized in the following equations. For simplicity we will limit the discussion to a single channel. The TOD, $`d`$ say, is modeled as the result of the action of a linear operator $`A`$ which encodes the optics and scanning strategy of the instrument, on the sky, made up of signal $`s`$ and foregrounds $`f`$:

$$d=A(s+f)+n.$$
(2)
Our assumptions is that $`s`$ is an isotropic Gaussian random field with zero mean and power spectrum $`C_{\mathrm{}}`$. We will take the noise correlations $`N=nn^T`$ and $`A`$ as given—though one of the statistical challenges of CMB analysis is to relax this assumption to some degree. We will not discuss this further here.

The task is then to extract as much information as possible about the cosmological parameters $`\mathrm{\Theta }`$ from the TOD. To set up this inverse problem we write down Bayes’ theorem

$$P(s,C_{\mathrm{}},f,\mathrm{\Theta }|d)P(d)=P(d|f,s,C_{\mathrm{}},\mathrm{\Theta })P(f,s,C_{\mathrm{}},\mathrm{\Theta }).$$
(3)
In the Bayesian context, solving the inverse problem means exploring and summarizing the posterior density $`P(s,C_{\mathrm{}},f,\mathrm{\Theta }|d)`$. On the right hand side we can use the conditional independence of the $`C_{\mathrm{}}`$ and the data given $`s`$, and the plausible independence of $`f`$ and $`s`$ to simplify

$$P(f,s,C_{\mathrm{}},\mathrm{\Theta })=P(f)P(s|C_{\mathrm{}})P(C_{\mathrm{}}|\mathrm{\Theta })P(\mathrm{\Theta }).$$
(4)
Traditionally, inference is performed in a linear sequence of steps which are concatenated into a pipeline. At each individual step a likelihood is written for the data in its current representation (e.g. TOD, maps, $`C_{\mathrm{}}`$) in terms of the parameters describing the next stage of compression. Due to the complexity of evaluating the likelihoods, often the likelihood approach is abandoned and approximate, suboptimal but unbiased estimators are constructed.

Within the above framework we can understand each of the steps as a limit of Eq. 3. For example the “map-making” step is signal estimation with

$$P(f,s,C_{\mathrm{}},\mathrm{\Theta })=P(f)P(s|C_{\mathrm{}})P(C_{\mathrm{}}|\mathrm{\Theta })P(\mathrm{\Theta })=const.$$
(5)
Similarly, “Power spectrum estimation” summarizes the marginal posterior $`P(C_{\mathrm{}}|d)`$, where the $`d`$ is taken in compressed form as the estimate of the CMB component. The summary proceeds through a maximum likelihood estimate (MLE) and the evaluation of the curvature of $`P(C_{\mathrm{}}|d)`$ at the MLE. In fact in most practical cases a non-optimal estimator is used due to the computational complexity of evaluating $`P(C_{\mathrm{}}|d)`$ or its derivatives and hence of computing the MLE.

Finally, “parameter estimation” signifies the step from the MLE of $`C_{\mathrm{}}`$ to a density $`P(\mathrm{\Theta }|d)`$. The posterior $`P(\mathrm{\Theta }|d)`$ is usually summarized in terms of the marginalized means and variances of the parameters or in terms of two-dimensional projections or marginalizations through the n-dimensional parameter space.

Recently, we developed MAGIC, a method that shows promise for global inference from the joint posterior MAGIC . The method, based on iterative sampling from the joint posterior, exploits the full dependence structure of the different science products from a CMB mission. For example it is useful for the separation of signal and foregrounds if the covariance of the signal is known. The information obtained on $`C_{\mathrm{}}`$ can be fed back into the CMB map estimate, which would in turn lead to a better estimate of $`C_{\mathrm{}}`$, etc. But since I discussed this technique in my contributed talk wandeltcontrib  I will not go into details here.

I will now discuss a selection of the challenges we face in the “map-making,” “power spectrum estimation” and “parameter estimation” steps on the way from the TOD to $`\mathrm{\Theta }`$.

## V Challenges in CMB Analysis

### V.1 Map-making Challenges

Pixelizations of the sphere. A very basic challenge in CMB analysis is the fact that the CMB is a random field on the sphere. Convenient numerical techniques for storing and manipulating functions on the sphere needed to be developed. For CMB analysis in particular, fast methods for spherical harmonic transforms had to be developed and implemented. Point sets were needed that had good generalized quadrature properties such that discrete numerical approximations to Eq. 1 and its inverse

$$a_\mathrm{}m=d^2\widehat{n}Y_\mathrm{}m(\widehat{n})T(\widehat{n})$$
(6)
generate accurate results. For local operations such as nearest neighbor searches and multi-resolution work, pixelizations of the sphere that allow hierarchical refinement are useful.

Various pixelizations of the sphere have been proposed HEALPix ; Igloo ; GLESP . Of these, HEALPix features pixels with exactly equal areas throughout, approximate equidistribution of pixel centers over the sphere, simple analytical equations for the pixel boundaries, fast spherical harmonic transforms and very favorable quadrature properties. Due to these features, HEALPix has developed into the standard pixelization for astrophysical all sky maps.

Beam Deconvolution Map-making. It is easy to show that the MLE $`\widehat{m}`$ for the map $`m=s+f`$ is the result of maximizing Eq. 3 using Eq. 5. This results in

$$m(A^TN^1A)^1A^TN^1d$$
(7)
with the associated noise covariance matrix $`C_N=<mm^T>=(A^TN^1A)^1`$.

These equations are generally valid whatever the forms of $`A`$ and $`N`$. However, the non-trivial structure in the observation matrix $`A`$ which is a consequence of the unavoidable imperfections in satellite optics makes solving for the map challenging. It is important to deconvolve the beam functions, since the beam distortions and side lobes lead to spurious shadow images of bright foregrounds. These images are spread throughout the map in a complicated way that depends on the scanning strategy. In addition, not accounting for beam imperfections results in distorting the signal itself. Both of these effects can bias not just the map but also the result of the covariance estimation.

For polarization map-making this deconvolution becomes even more important, since the polarization signal is weak and beam asymmetries can introduce spurious polarization signals into the data. Further, if bright foregrounds are significantly polarized, they may induce a significant polarization through their shadow images if beam convolution effects are neglected in map-making. Making high-quality maps of the polarization of the CMB anisotropy is the next frontier in CMB map-making.

To simulate the effects of realistic beams we have to have a general convolution technique for a beam map with a sky map along the scan path. If implemented as discrete sums over two pixelized maps, rotated to all possible relative orientations, such a technique requires of order $`n_p^{2.5}`$ operations, where $`n_p`$ is the number of pixels in the maps. By doing the convolution in spherical harmonic space it was shown in wandeltgorski  that this can be reduced to $`n_p^2`$ in the general case and to $`n_p^{3/2}`$ in interesting limit cases. This fast convolution method was generalized to polarization in challinor .

Component separation. Once a map of the sky is made from each channel, these maps can be compressed losslessly into maps at each frequency. A great deal of work in the field has gone into devising methods for then obtaining an estimate of the CMB sky from these foreground contaminated maps at each frequency band, both for temperature and for polarization. We can either choose to model the foregrounds physically (e.g. Stolyarov ; WMAPFG ) or we can attempt “blind separation,” by defining an algorithm for automatic detection of different components in the maps, e.g. based on the statistical independence of these components (e.g. Maino ).

Lensing. The CMB has traveled past cosmological mass concentrations which perturb the photon geodesics. This distortion contains valuable and complementary information about cosmological parameters. At the same time it mixes the polarization modes in the CMB data, contaminating the primary signal for the detection of the primordial gravitational wave background. Methods need to be devised that can measure this distortion and extract the information contained in it. It is intriguing that HirataSeljak  find that exact techniques have significant advantages over approximate, quadratic estimators for the reconstruction of B polarization maps from lensed CMB.

### V.2 Power Spectrum Estimation

The Computational Problem. For perfect (all-sky, pure signal, no noise) data, power spectrum estimation would be easy. Just compute the spherical harmonic transform, Eq. 6. Then the estimator

$$\widehat{C_{\mathrm{}}}=\frac{_m|a_{\mathrm{}}|^2}{2\mathrm{}+1}$$
(8)
is the MLE for $`C_{\mathrm{}}`$. It is also easy to evaluate and explore the perfect data posterior to quantify the uncertainty in the estimates. However, in the general case we would like to evaluate $`P(S(C_{\mathrm{}})|d,N)`$ which is the result of integrating out (“marginalizing over”) $`s`$ in the joint posterior. We obtain555We use $`G(x,X)`$ as a shorthand for the multivariate Gaussian density
$$G(x,X)=\frac{1}{\sqrt{|2\pi X|}}\mathrm{exp}\left(\frac{1}{2}x^TX^1x\right).$$
(9)

$$P(S(C_{\mathrm{}})|d,N)=G(m,S(C_{\mathrm{}})+C_N).$$
(10)
Here $`S`$ is the signal covariance matrix, parameterized by the $`C_{\mathrm{}}`$. Since $`S+C_N`$ is not a sparse matrix and since the determinant in the Gaussian depends on $`S`$, evaluating $`P(S(C_{\mathrm{}})|d,N)`$ as function of the $`C_{\mathrm{}}`$ costs of order $`n_p^3`$ operations.

For Planck, $`n_p10^7`$ so if the constant factor in the scaling law was $`1`$ (an unrealistic underestimate) this would mean $`10^{21}`$ operations for one likelihood evaluation. For a 10 GFLOP CPU this means 1000s of CPU years of computation.

Various approaches have been suggested to to counter this challenge. They can be broadly divided into two classes: 1) specialized exact (maximum likelihood estimation) algorithms exploit advantageous symmetries in the observational strategy of a CMB mission OSH ; WH  to reduce the computational scaling from order $`n_p^3`$ to order $`n_p^2`$, and 2) approximate algorithms which filter unwanted properties of the data and simply compute the power spectra on the filtered and incomplete data, and then de-bias the results using Monte Carlo simulations after the factmaster . This second class of algorithms is known as pseudo-$`C_{\mathrm{}}`$ algorithmsWHG  and has gained a great deal of popularity in recent analyses of CMB data, including the WMAP data WMAPCl .

Cosmic Variance. Aside from computational problems there is an interesting conceptual problem with CMB power spectrum analysis: the fact that we only have one sky.

This fact induces a fundamental limitation to how well we will be able to constrain cosmological parameter estimates from the CMB. The usual way of phrasing this limitation is in terms of cosmic variance.

Essentially, power spectrum estimation is variance estimation. For perfect data the solution that maximizes Eq. 10 is

$$\widehat{C_{\mathrm{}}}=\frac{_m|a_{\mathrm{}}|^2}{2\mathrm{}+1}.$$
(11)
So we are estimating the variance of the spherical harmonic coefficients. But, for example, for $`\mathrm{}=2`$ there are only 5 such coefficients and a variance estimate from 5 numbers is statistically uncertain. So “cosmic variance” is expression of the fundamental limit to the precision of any measurement of the $`C_{\mathrm{}}`$ caused by the fact that the sphere is a bounded space and the fact that causality will not allow us to observe independent patches of the universe.

This fundamental limit to our knowledge provides a powerful motivation to do the best possible job in analyzing cosmological data.

### V.3 Parameter Estimation

Techniques. The problem of parameter estimation is challenging since we need to explore the posterior density $`P(\mathrm{\Theta }|d)`$ which varies over $`1020`$ dimensions. Various techniques have been used to do this. For smaller number of dimensions (up to about 5) gridding techniques have worked well. However, the current state of the art is to use Markov Chain Monte Carlo methods Knox  such as the Metropolis Hastings algorithm to sample from $`P(\mathrm{\Theta }|d)`$ and to then base inferences on summaries of the posterior density computed from the sampled representation.

The question how to implement physical priors and constraints was one theme that was discussed at this conference. Within the Bayesian framework there exists a unique prescription for applying physical constraints through the specification of informative priors. The first example of a CMB parameter estimation which can be explored interactively online is the Cosmic Concordance Project HWC . This compiles data from several recent CMB observations and combines them with the user’s choice of other, non-CMB experiments and physical priors. The result is displayed as the 2D marginal posterior density for any 2 parameters chosen by the user. An example is displayed in figure 3. I invite you to have a look at our prototype implementation at http://galadriel.astro.uiuc.edu/ccp.

A first scientific result from this project was a measurement of the fraction of 4He in the Universe, both from CMB data alone and in the combination of standard Big Bang Nucleosynthesis with the CMB data. Since standard BBN links the primordial 4He abundance to the baryon to photon ratio which is determined exquisitely well by CMB data, we obtain the most precise measurement of the primordial 4He abundance to date HCW .

Several open questions remain to be addressed. What number $`n`$ of parameters is needed to fit the data? Which selection of parameters from the full set should we use in the analysis? These questions are currently handled on the basis of personal preference of the authors or computational convenience—a statistically motivated procedure for letting the data selecting certain parameters has not yet been implemented.

Operationally, it is very computationally expensive to sample from a posterior with 10 parameters or more, because each likelihood evaluation requires running a Boltzmann code which computes $`C(\mathrm{\Theta })_{\mathrm{}}`$. In many dimensions the Metropolis sampler produces correlated samples (regardless of whether the target density is correlated or not—the correlations come from the sampler taking small steps through the many dimensional space). New numerical techniques for evaluating and sampling from the likelihood in high-dimensional spaces are needed.

Beating cosmic variance. In Figure 2 we show another application of parameter estimation techniques which is relevant to several themes we touched on in this talk and at this conference. The figure shows our compilation of CMB data and the mean and $`\pm 68`$% errors on the range of power spectra which are contained in our 10 parameter fits. In other words, we have implemented the optimal non-linear filter for the $`C_{\mathrm{}}`$ if the Universe is really described by our 10 parameter model. The physical prior has reduced the cosmic variance error bars far below the limit set by Eq. 11. Since we are using the physical prior that the CMB power spectrum is the result of plasma oscillations in the primordial photon baryon fluid the resulting smoothness of the power spectrum is used in constructing the estimate.

## VI Testing the Assumptions: Challenges of Non-canonical CMB Analysis

In addition to the canonical analyses outlined above it is fundamentally important to test the assumptions on which these analyses rest. I will very briefly mention two of them here.

Is the CMB signal really Gaussian? Is it isotropic? These questions touch simultaneously on the issue of hypothesis testing, as well as model selection. What is the evidence in the data for the assumptions of isotropy and Gaussianity on which canonical analyses of the CMB are built? Can a goodness of fit criterion be defined which allows assessing whether the standard model of cosmology is a complete description of the CMB data?

The idea of testing the goodness of fit can be extended to cross tests of CMB data with other cosmological data. Ultimately the goal of the Cosmic Concordance Project is to allow users to select various data sets and explore interactively whether the parameter constraints from various data sets are compatible with each other. If this agreement is established the the constraints can be combined to generate yet stronger constraints. If disagreement is found this is motivation for observational groups to collect more data or for theoretical groups to work out new mechanisms that can reconcile the discrepant observations.

Testing for statistical isotropy in the CMB is a well-defined operation, since it is easy to specify alternative models. It requires checking whether a model of the correlations in the fluctuations in terms of rotationally invariant quantities (such as the power spectrum) is better or worse than a model that contains quantities that do not transform as scalars under SO(3). A frequentist test statistic has been suggested in hajiansouradeep , but a Bayesian treatment has not yet been attempted. The detection of a significant deviation from statistical isotropy would be a very important result, since isotropy is a fundamental prediction of inflation.

Testing for non-Gaussianity (NG) is similarly important, but much less well-defined. In the absence of physical models for NG a Bayesian treatment is not possible. The standard approach is to define some NG statistic, e.g. the skewness of the one-point function of the CMB fluctuations. Then this statistic is applied to the data and to a sample of Gaussian Monte Carlo samples (which are usually constrained to match the two-point statistics of the observed data). This Monte Carlo sample represents the null-hypothesis. A way is defined to assess discrepancy of the data and the Monte Carlo sample and if the discrepancy is statistically significant a detection is claimed.

Even though this is a straightforward frequentist procedure there is a great deal of arbitrariness in choosing the test statistic. Usually the choice is made based on a vague notion of genericity, in the sense that, for example, the skewness of the data is probably a more generic NG statistic than the temperature in pixel number 2,437,549, say. However, it is clearly not trivial to define a proper measure on the space of all possible statistics. The arbitrariness in the selection of the statistic (topological features of the map, n-point functions in pixel and in spherical harmonic space, wavelets, excursion sets, etc.) makes the test results difficult to interpret. For example, is clearly true that there will always be some statistic that gives an n-sigma result for any n!

To circumvent this arbitrariness, we could test for the NG predicted in certain variants of inflationary models (KomatsuSpergel ; KSW ; WMAPNG  and references therein). These physical models allow a Bayesian analysis in principle, and it would provide the best possible constraints on the presence of NG, but implementing Bayesian NG inference these models is computationally tedious.

## VII Conclusion

In this talk I reviewed some of the more severe statistical and computational challenges of CMB analysis. CMB data are fundamentally important to cosmology, the problems I outlined are intellectually rich, and the knowledge that we need to make the most of the one CMB sky we can observe and analyze, present a powerful motivation to solve the CMB analysis problem.

###### Acknowledgements.

I thank the organizers of PHYSTAT2003 for making this conference so enjoyable and stimulating. This work has been partially supported by the National Computational Science Alliance under grant number AST020003N and by the University of Illinois at Urbana-Champaign through an NCSA Faculty Fellowship.