text:
# Möbius Symmetry of Discrete Time Soliton Equations

## 1 Introduction

It has been known that, for a given soliton equation, there exists a discrete analogue of the evolution equation, which preserves integrability. Once we find a discrete integrable equation, we can derive infinitely many integrable differential equations by taking continuous limits of variables in many different ways. We are interested in characterizing such integrable discrete systems.

In our previous papers  we studied various integrable discrete systems which have certain symmetry and satisfy periodic boundary conditions. In this approach the nature of a quadratic equation, imposed by the boundary condition, plays a crucial role. The map which generates a time evolution turns out to be twofold, corresponding to two solutions of the quadratic equation. We showed that the discriminant of the quadratic equation becomes a perfect square when the system is integrable. Owing to this fact the map is free from a square root, hence is rational. We will call this type of map a non square root map (NSRM) in what follows. The discrete versions of Lotka-Volterra equation, KdV equation, Toda lattice, KP equation and Painlevé equations  are such examples. The purpose of this paper is to generalize this scheme of characterizing discrete integrable systems.

To clarify the point of our argument let us show first a typical example. The 3-point discrete time Toda lattice is given by the set of coupled equations:

$$\{\begin{array}{ccc}\overline{x}_k\overline{u}_k& =& u_kx_{k1}\\ \overline{x}_k+\overline{u}_{k+1}& =& u_k+x_k\end{array},k=1,2,3.$$
(1)
Here $`\overline{x}_k`$ means the variable $`x_k`$ at the time $`t+1`$, i.e., $`\overline{x}_k(t)=x_k(t+1)`$. We impose the periodic conditions $`x_{k+3}=x_k`$ and $`u_{k+3}=u_k`$. Hence there are 6 dependent variables.

Before solving eq.(1) directly it is more instructive to solve

$$\{\begin{array}{ccc}\overline{x}_k\overline{u}_k& =& p_k\\ \overline{x}_k+\overline{u}_{k+1}& =& q_k\end{array}$$
(2)
for arbitrary functions $`p_k,q_k`$ of $`(x,u)`$. We can also write them as

$$\{\begin{array}{ccc}\overline{x}_k& =& \frac{q_k\overline{x}_{k+1}p_{k+1}}{\overline{x}_{k+1}}\\ \overline{u}_k& =& \frac{q_{k1}\overline{u}_{k1}p_{k1}}{\overline{u}_{k1}}.\end{array}$$
(3)
Note that they are special cases of Möbius transformation. Hence by repeating this transformation, say for the first equation of eq.(3), three times, $`\overline{x}_k`$ is given in terms of $`(\overline{x}_{k+3},\overline{u}_{k3})`$, which are nothing but $`(\overline{x}_k,\overline{u}_k)`$ themselves owing to the periodic boundary conditions. In this way we obtain an equation

$$\overline{x}_k=\frac{A\overline{x}_k+B}{\mathrm{\Gamma }\overline{x}_k+\mathrm{\Delta }}$$
(4)
or, equivalently,

$$\mathrm{\Gamma }\overline{x}_k^2(A\mathrm{\Delta })\overline{x}_kB=0,$$
(5)
where

$`A`$ $`=`$ $`q_1q_2q_3q_kp_kq_{k1}p_{k1},`$
$`B`$ $`=`$ $`q_kq_{k+1}p_{k+1}+p_{k+1}p_{k1},`$
$`\mathrm{\Gamma }`$ $`=`$ $`q_{k+1}q_{k1}p_k,`$
$`\mathrm{\Delta }`$ $`=`$ $`q_{k+1}p_{k+1},`$
in the case of three point Toda lattice.

The quadratic equation eq.(5) has two solutions

$$\overline{x}_k=\frac{A\mathrm{\Delta }\pm \sqrt{Dis}}{2\mathrm{\Gamma }},$$
(6)
unless its discriminant

$$Dis=(A\mathrm{\Delta })^2+4B\mathrm{\Gamma }$$
(7)
vanishes. The general form of solutions contain square root of polynomial functions. Under this circumstance a sequence of the map will yield very complicated orbits in general. But the integrable systems are not such cases, i.e., the discriminant happens to be a perfect square of some polynomial. In fact if we substitute the right hand sides of eq.(1) into $`p_k`$’s and $`q_k`$’s of the expression eq.(1), the discriminant is considerably simplified and is given by

$$Dis=(x_1x_2x_3u_1u_2u_3)^2.$$
(8)
This is a typical non square root map (NSRM), which we are going to study (see below for more precise definition).

This remarkable feature could be understood from other point of view as well. Namely by an inspection of our original equations eq.(1) we can see easily that

$$\overline{x}_k=u_k,\overline{u}_k=x_{k1}$$
(9)
is one of sets of solution. This implies that another set of solutions of the quadratic equation eq.(5) must not have a square root either. We would like to emphasize that the existence of the solution of such simple form eq.(9) owes to the symmetry of the equation eq.(1).

We will study, in the following section, $`N`$-point maps under which the square roots are removed when periodic boundary conditions are imposed. In section 3 the discrete time Toda equation is studied in detail based on the preceding argument. In particular we will show that this scheme enables one to generalize the Toda equation to a $`q`$-difference form, so that various $`q`$-difference soliton equations are derived by reductions, including the $`qP_{IV}`$ equation by Kajiwara, Noumi and Yamada. In section 4, we will discuss time evolution of solutions. The last section is devoted to summarize arguments in this paper.

## 2 $`N`$-point Map with Periodic Boundary Conditions

We consider $`N`$-point systems in this section and study conditions under which a map becomes a non square root map (NSRM).

### 2.1 Möbius transformation

Generalizing eq.(3), let us consider the following sequence of Möbius transformations,

$$\overline{x}_{k+h(n1)}=\frac{\alpha _n\overline{x}_{k+hn}+\beta _n}{\gamma _n\overline{x}_{k+hn}+\delta _n},n=1,2,\mathrm{},N,h=\pm 1,\alpha _n\delta _n\beta _n\gamma _n0$$
(10)
and Möbius matrices

$$M_n=\left(\begin{array}{cc}\alpha _n& \beta _n\\ \gamma _n& \delta _n\end{array}\right),detM_n0,n=1,2,\mathrm{},N.$$
(11)
After $`N`$ times of the transformation we obtain the equation

$$\overline{x}_k=\frac{A\overline{x}_{k+hN}+B}{\mathrm{\Gamma }\overline{x}_{k+hN}+\mathrm{\Delta }},$$
(12)
while the corresponding Möbius map is given by

$$S_N:=\left(\begin{array}{cc}A& B\\ \mathrm{\Gamma }& \mathrm{\Delta }\end{array}\right)=M_1M_2\mathrm{}M_N.$$
(13)
We now impose the periodic boundary conditions $`\overline{x}_k=\overline{x}_{k+hN}`$. Then the formula eq.(12) is equal to the quadratic equation

$$\mathrm{\Gamma }\overline{x}_k^2(A\mathrm{\Delta })\overline{x}_kB=0.$$
(14)
The solutions of this equation are given by

$$\overline{x}_k=\frac{\mathrm{tr}S_N2\mathrm{\Delta }\pm \sqrt{Dis}}{2\mathrm{\Gamma }}$$
(15)
whereas the discriminant $`Dis`$ is

$$Dis=(A\mathrm{\Delta })^2+4B\mathrm{\Gamma }=(\mathrm{tr}S_N)^24detS_N.$$
(16)
### 2.2 Non square root maps

In this subsection we study maps such that the discriminant eq.(16) becomes a perfect square of some polynomial function. As we will prove shortly the following 6 patterns of the map $`M_n`$ satisfy the requirement:

$$\left(\begin{array}{cc}0& \beta _n\\ \beta _n& 0\end{array}\right),\left(\begin{array}{cc}\alpha _n& 0\\ 0& \delta _n\end{array}\right),\left(\begin{array}{cc}\alpha _n& \beta _n\\ 0& \delta _n\end{array}\right),$$

$$\left(\begin{array}{cc}\alpha _n& 0\\ \gamma _n& \delta _n\end{array}\right),\left(\begin{array}{cc}\gamma _n+\beta _n& \beta _n\\ \gamma _n& 0\end{array}\right),\left(\begin{array}{cc}0& \beta _n\\ \gamma _n& \beta _n+\gamma _n\end{array}\right).$$
(17)
Proof:

Since the trace and the determinant are invariant under the transposition ot the matrices, we do not consider here the fourth and sixth patterns. For the first three patterns we find

$$S_N=\beta \left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)^N,\left(\begin{array}{cc}\alpha & 0\\ 0& \delta \end{array}\right),\left(\begin{array}{cc}\alpha & F\\ 0& \delta \end{array}\right)$$
(18)
where

$$\alpha =\alpha _1\mathrm{}\alpha _N,\beta =\beta _1\mathrm{}\beta _N,\delta =\delta _1\mathrm{}\delta _N,$$

and $`F`$ is a certain function of $`\alpha _n,\beta _n,\delta _n`$. From this expression we easily find

$$Dis=(2\beta )^2\delta _{N,odd},(\alpha \delta )^2,(\alpha \delta )^2$$

corresponding to the first, second and the third patterns.

For the 5th pattern we first observe

$$detS_N=detM_1\mathrm{}detM_N=\beta \gamma .$$
(19)
In order to calculate $`\mathrm{tr}S_N`$, we split $`M_n`$ into two parts

$$M_n=\left(\begin{array}{cc}\beta _n& \beta _n\\ 0& 0\end{array}\right)+\left(\begin{array}{cc}\gamma _n& 0\\ \gamma _n& 0\end{array}\right)B_n+C_n.$$
(20)
Using the property

$$B_mC_n=0,m,n1,2,\mathrm{}N,$$

we can see that $`S_N`$ has terms of the form $`C_1\mathrm{}C_\rho B_{\rho +1}\mathrm{}B_N`$ only. From this fact, the trace of $`S_N`$ is calculated as

$`\mathrm{tr}S_N`$ $`=`$ $`{\displaystyle \underset{\rho N}{}}\mathrm{tr}(C_1\mathrm{}C_\rho B_{\rho +1}\mathrm{}B_N)`$
$`=`$ $`{\displaystyle \underset{\rho N}{}}\mathrm{tr}(C_2\mathrm{}C_\rho B_{\rho +1}\mathrm{}B_NC_1)`$
$`=`$ $`\mathrm{tr}(B_1\mathrm{}B_N)+\mathrm{tr}(C_1\mathrm{}C_N)`$
$`=`$ $`\beta +\gamma `$
where $`\gamma =\gamma _1\mathrm{}\gamma _N`$. Similarly we can treat the 6th pattern and obtain the same results as eqs.(19) and (2.2) for the determinant and trace, respectively.

Hence we have obtained discriminants for all six patterns as follows:

| 1) | $`\mathrm{tr}S_N=\delta _{N,even}\beta `$, | $`detS_N=\beta ^2`$ | $``$ | $`Dis=(2\beta )^2\delta _{N,odd}`$, |
| --- | --- | --- | --- | --- |
| 2) | $`\mathrm{tr}S_N=\alpha +\delta `$, | $`detS_N=\alpha \delta `$ | $``$ | $`Dis=(\alpha \delta )^2`$, |
| 3) | $`\mathrm{tr}S_N=\alpha +\delta `$, | $`detS_N=\alpha \delta `$ | $``$ | $`Dis=(\alpha \delta )^2`$, |
| 4) | $`\mathrm{tr}S_N=\alpha +\delta `$, | $`detS_N=\alpha \delta `$ | $``$ | $`Dis=(\alpha \delta )^2`$, |
| 5) | $`\mathrm{tr}S_N=\gamma +\beta `$, | $`detS_N=\beta \gamma `$ | $``$ | $`Dis=(\gamma \beta )^2`$, |
| 6) | $`\mathrm{tr}S_N=\gamma +\beta `$, | $`detS_N=\beta \gamma `$ | $``$ | $`Dis=(\gamma \beta )^2`$. |

Therefore every discrete equation with the Möbius matrix pattern of eq.(17) has no square root, although they are not all.

### 2.3 Rational maps

Solutions of the quadratic equation eq.(14) are given by eq.(15). The purpose of this subsection is to write them explicitly in the case of NSRM. First we note that the 2nd and the 3rd pattern as well as even $`N`$ case of the 1st pattern do not form a quadratic equation due to $`\mathrm{\Gamma }=0`$. Hence it is enough to consider the odd $`N`$ case of 1st pattern, 4th, 5th and 6th patterns.

#### 2.3.1 Odd $`N`$ 1st pattern

The odd $`N`$ case is given by the following setup:

$$S_N=\left(\begin{array}{cc}0& \beta \\ \beta & 0\end{array}\right),\mathrm{tr}S_N=0,detS_N=\beta ^2,Dis=(2\beta )^2.$$
(22)
Hence the solutions of the quadratic equation are

$$\overline{x}_k=\frac{\mathrm{tr}S_N2\mathrm{\Delta }\pm \sqrt{Dis}}{2\mathrm{\Gamma }}=\frac{\pm 2\beta }{2\beta }=\pm 1.$$
(23)
#### 2.3.2 4th pattern

In this pattern we have

$$S_N=\left(\begin{array}{cc}\alpha & 0\\ F& \delta \end{array}\right),\mathrm{tr}S_N=\alpha +\delta ,detS_N=\alpha \delta ,Dis=(\alpha \delta )^2.$$
(24)
To calculate $`F`$, we split $`M_n`$ to two parts $`A_n`$ and $`D_n`$ as follows:

$$M_n=\left(\begin{array}{cc}\alpha _n& 0\\ \gamma _n& \delta _n\end{array}\right)=\left(\begin{array}{cc}\alpha _n& 0\\ 0& 0\end{array}\right)+\left(\begin{array}{cc}0& 0\\ \gamma _n& \delta _n\end{array}\right)A_n+D_n.$$
(25)
Using the fact $`A_mD_n=0`$, we see that $`S_N`$ has terms of the form $`D_1\mathrm{}D_{\rho 1}A_\rho \mathrm{}A_N`$ only, where $`\rho =1,2,\mathrm{},N+1`$ and $`D_0=A_{N+1}=0`$. All possible terms can be written in such forms as

$`A_1A_2\mathrm{}A_N`$ $`=`$ $`\left(\begin{array}{cc}\alpha & 0\\ 0& 0\end{array}\right),`$ (28)
$`D_1D_2\mathrm{}D_N`$ $`=`$ $`\left(\begin{array}{cc}0& 0\\ \delta _1\mathrm{}\delta _{N1}\gamma _N& \delta \end{array}\right),`$ (31)
$`D_1\mathrm{}D_{\rho 1}A_\rho \mathrm{}A_N`$ $`=`$ $`\left(\begin{array}{cc}0& 0\\ \delta _1\mathrm{}\delta _{\rho 2}\gamma _{\rho 1}\alpha _\rho \mathrm{}\alpha _N& 0\end{array}\right).`$ (34)
Hence we find the expression of $`F`$ as

$$F=\underset{\rho =2}{\overset{N}{}}\delta _1\mathrm{}\delta _{\rho 2}\gamma _{\rho 1}\alpha _\rho \mathrm{}\alpha _N.$$
(35)
The solutions of the quadratic equation are then given as follows:

$$\overline{x}_k=\frac{\mathrm{tr}S_N2\mathrm{\Delta }\pm \sqrt{Dis}}{2\mathrm{\Gamma }}=\frac{(\alpha \delta )\pm (\alpha \delta )}{2F}=\{\begin{array}{c}\frac{\alpha \delta }{F}\\ \\ 0\end{array}.$$
(36)
#### 2.3.3 5th pattern

We already got $`\mathrm{tr}S_N`$ and $`detS_N`$ in the previous subsection, but did not have calculated the elements of $`S_N`$. To obtain the solution of the quadratic equation, we must find their elements since information about $`C`$ and $`D`$ are necessary for the solutions.

The splitting idea of eq.(20) is convenient to obtain the elements here again. Using this idea , we need to calculate only the following quantities

$`C_1C_2\mathrm{}C_N`$ $`=`$ $`\gamma \left(\begin{array}{cc}1& 0\\ 1& 0\end{array}\right),`$ (39)
$`B_1B_2\mathrm{}B_N`$ $`=`$ $`\beta \left(\begin{array}{cc}1& 1\\ 0& 0\end{array}\right),`$ (42)
$`C_1\mathrm{}C_{\rho 1}B_\rho \mathrm{}B_N`$ $`=`$ $`\gamma _1\mathrm{}\gamma _{\rho 1}\beta _\rho \mathrm{}\beta _N\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),`$ (45)
because $`B_mC_n=0`$. Hence $`\mathrm{\Gamma }`$ and $`\mathrm{\Delta }`$ of the elements of $`S_N`$ are obtained as

$`\mathrm{\Gamma }`$ $`=`$ $`\gamma _1\mathrm{}\gamma _N+{\displaystyle \underset{\rho =2}{\overset{N}{}}}\gamma _1\gamma _2\mathrm{}\gamma _{\rho 1}\beta _\rho \mathrm{}\beta _N,`$ (46)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \underset{\rho =2}{\overset{N}{}}}\gamma _1\gamma _2\mathrm{}\gamma _{\rho 1}\beta _\rho \mathrm{}\beta _N.`$ (47)
Consequently, noting the fact $`\mathrm{\Gamma }=\gamma \mathrm{\Delta }`$, we obtain the solutions of the quadratic equation as

$`\overline{x}_k`$ $`=`$ $`{\displaystyle \frac{\mathrm{tr}S_N2\mathrm{\Delta }\pm \sqrt{Dis}}{2\mathrm{\Gamma }}}={\displaystyle \frac{(\gamma +\beta )2\mathrm{\Delta }\pm (\gamma \beta )}{2(\gamma \mathrm{\Delta })}}`$ (51)
$`=`$ $`\{\begin{array}{ccc}1& & \\ & & \\ {\displaystyle \frac{\beta \mathrm{\Delta }}{\gamma \mathrm{\Delta }}}& & {\displaystyle \frac{\beta _N}{\gamma _1}}{\displaystyle \frac{Q_N(\beta ,\gamma )}{Q_1(\beta ,\gamma )}}\end{array}`$
where $`Q_n(\beta ,\gamma )`$ is the homogeneous polynomial of $`(\beta _1,\mathrm{},\beta _N,\gamma _1,\mathrm{},\gamma _N)`$, that is

$$Q_n(\beta ,\gamma )=Q_n(\beta _1,\mathrm{},\beta _N,\gamma _1,\mathrm{},\gamma _N)=\underset{\rho =2}{\overset{N}{}}\gamma _{n+1}\mathrm{}\gamma _{n+\rho 1}\beta _{n+\rho }\mathrm{}\beta _{n+N1}.$$
(52)
#### 2.3.4 6th pattern

For this pattern too, we already got $`\mathrm{tr}S_N`$ and $`detS_N`$ in the previous subsection, but did not have calculated the elements of $`S_N`$. Using the splitting idea again, the elements of $`S_N`$ in this pattern are able to be obtained. The splitting of $`M_n`$ in this case is

$$M_n=\left(\begin{array}{cc}0& \beta _n\\ \gamma _n& \beta _n+\gamma _n\end{array}\right)=\left(\begin{array}{cc}0& \beta _n\\ 0& \beta _n\end{array}\right)+\left(\begin{array}{cc}0& 0\\ \gamma _n& \gamma _n\end{array}\right)B_n+C_n.$$
(53)
Since $`C_mB_n=0`$, we are enough to calculate the following terms

$`C_1C_2\mathrm{}C_N`$ $`=`$ $`\gamma \left(\begin{array}{cc}0& 0\\ 1& 1\end{array}\right),`$ (56)
$`B_1B_2\mathrm{}B_N`$ $`=`$ $`\beta \left(\begin{array}{cc}0& 1\\ 0& 1\end{array}\right),`$ (59)
$`B_1\mathrm{}B_{\rho 1}C_\rho \mathrm{}C_N`$ $`=`$ $`\beta _1\beta _2\mathrm{}\beta _{\rho 1}\gamma _\rho \mathrm{}\gamma _N\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right).`$ (62)
Hence we obtain $`\mathrm{\Gamma }`$ and $`\mathrm{\Delta }`$ of the elements of $`S_N`$ as

$`\mathrm{\Gamma }`$ $`=`$ $`\gamma {\displaystyle \underset{\rho =2}{\overset{N}{}}}\beta _1\mathrm{}\beta _{\rho 1}\gamma _\rho \mathrm{}\gamma _N,`$ (63)
$`\mathrm{\Delta }`$ $`=`$ $`\beta +\gamma +{\displaystyle \underset{\rho =2}{\overset{N}{}}}\beta _1\mathrm{}\beta _{\rho 1}\gamma _\rho \mathrm{}\gamma _N.`$ (64)
If we further use $`\mathrm{\Gamma }=\beta \mathrm{\Delta }`$, we get solutions of the quadratic equation as follows:

$`\overline{x}_k`$ $`=`$ $`{\displaystyle \frac{\mathrm{tr}S_N2\mathrm{\Delta }\pm \sqrt{Dis}}{2\mathrm{\Gamma }}}={\displaystyle \frac{(\gamma +\beta )2\mathrm{\Delta }\pm (\beta \gamma )}{2(\beta \mathrm{\Delta })}}`$ (68)
$`=`$ $`\{\begin{array}{ccc}1& & \\ & & \\ {\displaystyle \frac{\mathrm{\Delta }\gamma }{\mathrm{\Delta }\beta }}& & {\displaystyle \frac{\beta _1}{\gamma _N}}{\displaystyle \frac{Q_1(\gamma ,\beta )}{Q_N(\gamma ,\beta )}}.\end{array}`$
### 2.4 Polynomial $`Q`$

We note that polynomial $`Q`$’s, which appear in the solutions of the 5th and 6th patterns, are given by the minor determinants of the diagonal entries of the following matrix.

$$Q(\beta _1,\mathrm{},\beta _N,\gamma _1,\mathrm{},\gamma _N):=\left(\begin{array}{ccccccc}\beta _1+\gamma _1& \beta _1& 0& 0& \mathrm{}& 0& \gamma _1\\ \gamma _2& \beta _2+\gamma _2& \beta _2& 0& \mathrm{}& 0& 0\\ 0& \gamma _3& \beta _3+\gamma _3& & & & \mathrm{}\\ \mathrm{}& & & & & & 0\\ 0& 0& \mathrm{}& & & & \beta _{N1}\\ \beta _N& 0& \mathrm{}& & 0& \gamma _N& \beta _N+\gamma _N\end{array}\right).$$
(69)
For simplicity, we will introduce the following notations.

$$\begin{array}{ccc}Q(\beta ,\gamma )\hfill & :=& Q(\beta _1,\mathrm{},\beta _N,\gamma _1,\mathrm{},\gamma _N),\hfill \\ Q(\gamma ,\beta )\hfill & :=& Q(\gamma _1,\mathrm{},\gamma _N,\beta _1,\mathrm{},\beta _N),\hfill \\ Q^{}(\beta ,\gamma )\hfill & :=& Q(\beta _N,\mathrm{},\beta _1,\gamma _N,\mathrm{},\gamma _1),\hfill \\ Q^{}(\gamma ,\beta )\hfill & :=& Q(\gamma _N,\mathrm{},\gamma _1,\beta _N,\mathrm{},\beta _1).\hfill \end{array}$$

Namely, the matrix $`Q`$ having entries in ascending order is expressed by $`Q`$, and in descending order by $`Q^{}`$. These matrices are related by

$$Q^{}(\beta ,\gamma )=JQ(\gamma ,\beta )J,Q(\beta ,\gamma )=JQ^{}(\gamma ,\beta )J$$

where

$$J:=\left(\begin{array}{cccccccccc}0& \mathrm{}& \mathrm{}& 0& 1& & & & & \\ \mathrm{}& & & 1& 0& & & & & \\ \mathrm{}& & \mathrm{}& & \mathrm{}& & & & & \\ 0& 1& & & \mathrm{}& & & & & \\ 1& 0& \mathrm{}& \mathrm{}& 0& & & & & \end{array}\right).$$
(70)
We set symbols of them of $`(n,n)`$-entry as $`Q_n`$, i.e.

$$\begin{array}{cc}Q_n(\beta ,\gamma )& =\underset{\rho =1}{\overset{N}{}}\underset{\mu =1}{\overset{\rho 1}{}}\gamma _{n+\mu }\underset{\nu =\rho }{\overset{N1}{}}\beta _{n+\nu },\hfill \\ Q_n(\gamma ,\beta )& =\underset{\rho =1}{\overset{N}{}}\underset{\mu =1}{\overset{\rho 1}{}}\beta _{n+\mu }\underset{\nu =\rho }{\overset{N1}{}}\gamma _{n+\nu },\hfill \\ Q_n^{}(\beta ,\gamma )& =\underset{\rho =1}{\overset{N}{}}\underset{\mu =1}{\overset{\rho 1}{}}\gamma _{n+N\mu }\underset{\nu =\rho }{\overset{N1}{}}\beta _{n+N\nu },\hfill \\ Q_n^{}(\gamma ,\beta )& =\underset{\rho =1}{\overset{N}{}}\underset{\mu =1}{\overset{\rho 1}{}}\beta _{n+N\mu }\underset{\nu =\rho }{\overset{N1}{}}\gamma _{n+N\nu },\hfill \end{array}$$

where $`\beta _{n\pm N}=\beta _n,\gamma _{n\pm N}=\gamma _n`$.

These polynomials satisfy the following useful relations.

1. $`Q_{Nn+1}(\beta ,\gamma )=Q_n^{}(\gamma ,\beta ),Q_n(\beta ,\gamma )=Q_{Nn+1}^{}(\gamma ,\beta )`$,
2. $`(\beta _n+\gamma _n)Q_n(\beta ,\gamma )=\beta _{n1}Q_{n1}(\beta ,\gamma )+\gamma _{n+1}Q_{n+1}(\beta ,\gamma )`$.

The relation 1 is obvious. The relation 2 is non-trivial and will be checked by means of more fundamental relation :

$$\beta _nQ_n(\beta ,\gamma )=\beta \gamma +\gamma _{n+1}Q_{n+1}(\beta ,\gamma )$$

where $`\beta =\beta _1\mathrm{}\beta _N,\gamma =\gamma _1\mathrm{}\gamma _N`$.

We can show that $`Q_n(\beta ,\gamma )`$ coincides with the cofactor of the $`(n,n)`$ element of the product $`RL`$ of the following matrices

$$L=\left(\begin{array}{ccccccccc}1& 0& 0& \mathrm{}& \beta _N& & & & \\ \beta _1& 1& & 0& & & & & \\ 0& \beta _2& 1& \mathrm{}& & & & & \\ & & & \mathrm{}& & & & & \\ 0& 0& \mathrm{}& \beta _{N1}& 1& & & & \end{array}\right),R=\left(\begin{array}{ccccccccc}\gamma _1& 1& 0& \mathrm{}& 0& & & & \\ 0& \gamma _2& 1& & & & & & \\ & & \mathrm{}& & & & & & \\ 0& & & & 1& & & & \\ 1& 0& \mathrm{}& 0& \gamma _N& & & & \end{array}\right)$$

which could represent the dToda equations .

### 2.5 Modification

The discriminant $`Dis`$ is invariant under arbitrary similarity transformations, since it is given by the trace and determinant of $`S_N`$. This fact enables us to modify our map starting from the six patterns we considered above.

Namely under the transformation of $`M_n`$,

$$M_nM_n^{mod}U_n^1M_nU_{n+1},$$
(71)
the matrix eq.(13) is changed into

$$S_N^{mod}=M_1^{mod}M_2^{mod}\mathrm{}M_N^{mod}=U_1^1M_1M_2\mathrm{}M_NU_{N+1}.$$
(72)
If we take into account the periodic boundary condition, we must have $`U_{N+1}=U_1`$. Therefore the trace and the determinant of $`S_N`$, hence the discriminant $`Dis`$ as well, remain constant:

$$\mathrm{tr}S_N^{mod}=\mathrm{tr}S_N,detS_N^{mod}=detS_N,Dis^{mod}=Dis,$$
(73)
under the transformation eq.(71) for arbitrary matrices $`U_n`$’s. This implies that the modified map $`M_n^{mod}`$ is NSRM if $`M_n`$ is so.

Now denoting by $`\overline{x}^{mod}`$ the modified variable of $`\overline{x}`$, we will write the Möbius map associated with $`M_n`$ as $`M_n(x_k)`$. In this way, the modified Möbius map is possible to be written as

$$\overline{x}_k^{mod}=M_1^{mod}(\overline{x}_{k+1}^{mod})=U_1^1M_1U_1(\overline{x}_{k+1}^{mod}).$$
(74)
Here if we define $`\overline{x}_{k+n}^{mod}=U_{n+1}^1(\overline{x}_{k+n})`$, then

$$\begin{array}{cccc}& \overline{x}_k^{mod}& =& U_1^1M_1U_2(\overline{x}_{k+1}^{mod})\hfill \\ & U_1^1(\overline{x}_k)& =& U_1^1M_1U_2U_2^1(\overline{x}_{k+1})\hfill \\ & \overline{x}_k& =& M_1(\overline{x}_{k+1}).\hfill \end{array}$$
(75)
As a result, the following two modifications are equivalent

$$\begin{array}{cccccc}1)& M_n\hfill & & M_n^{mod}\hfill & =& U_n^1M_nU_{n+1},\hfill \\ 2)& \overline{x}_{k+n}\hfill & & \overline{x}_{k+n}^{mod}\hfill & =& U_{n+1}^1(\overline{x}_{k+n}).\hfill \end{array}$$
(76)
This correspondence is preserved in the solutions of the quadratic equation, so we can create many equations from the basic six patterns in this way.

## 3 Discrete Time Toda Equation

We study, in this section, the $`N`$-point discrete time Toda equation (dToda) in detail from the view point of the previous section.

### 3.1 Two types of solution

Let us start with writing the $`N`$-point dToda equation.

$$\{\begin{array}{ccc}\overline{x}_k\overline{u}_k& =& u_kx_{k1}\\ \overline{x}_k+\overline{u}_{k+1}& =& u_k+x_k\end{array},k=1,2,3,\mathrm{},N.$$
(77)
In order to apply the argument in the previous section to this system we transform eq.(77) into the form of Möbius transformation,

$`\{\begin{array}{ccc}\overline{x}_k\hfill & =& {\displaystyle \frac{(x_k+u_k)\overline{x}_{k+1}x_ku_{k+1}}{\overline{x}_{k+1}}}\hfill \\ \overline{u}_k\hfill & =& {\displaystyle \frac{(x_{k1}+u_{k1})\overline{u}_{k1}x_{k2}u_{k1}}{\overline{u}_{k1}}}\hfill \end{array}`$ (80)
corresponding to $`(\overline{x}_{k+1},\overline{u}_{k1})(\overline{x}_k,\overline{u}_k)`$.

We notice that this is a modified 5th pattern Möbius transformation which is specified by the following data:

$$(\beta _n,\gamma _n)=\{\begin{array}{c}(x_{k+n1},u_{k+n1})\\ (u_{kn},x_{kn})\end{array},U_n=\left(\begin{array}{cc}1& 0\\ 0& \gamma _n\end{array}\right).$$
(81)
Substituting these data into eq.(51) we obtain solutions immediately. Before writing them down, however, it is convenient to simplify notations. Namely the matrix $`Q`$, which appears in the solutions, is $`Q(x_k,x_{k+1},\mathrm{},x_{k2},x_{k1},u_k,u_{k+1},\mathrm{},u_{k2},u_{k1})`$. Let the matrix $`Q`$ be

$$\begin{array}{ccc}Q(x_k,u_k)\hfill & :=\hfill & Q(x_k,\mathrm{},x_{k+N1},u_k,\mathrm{},u_{k+N1}),\hfill \\ Q^{}(x_k,u_k)\hfill & :=\hfill & Q(x_{k+N1},\mathrm{},x_k,u_{k+N1},\mathrm{},u_k),\hfill \end{array}$$

and its minor determinants be

$$\begin{array}{ccc}Q_n(x_k,u_k)\hfill & :=\hfill & Q_n(x_k,\mathrm{},x_{k+N1},u_k,\mathrm{},u_{k+N1}),\hfill \\ Q_n^{}(x_k,u_k)\hfill & :=\hfill & Q_n(x_{k+N1},\mathrm{},x_k,u_{k+N1},\mathrm{},u_k).\hfill \end{array}$$
(82)
These minor determinants satisfy

1. $`Q_{Nn+1}(x_k,u_k)=Q_n^{}(u_k,x_k),Q_n(x_k,u_k)=Q_{Nn+1}^{}(u_k,x_k)`$,
2. $`(x_{k+n}+u_{k+n})Q_n(x_k,u_k)=x_{k+n1}Q_{n1}(x_k,u_k)+u_{k+n+1}Q_{n+1}(x_k,u_k),`$
3. $`Q_n(x_{k+N},u_{k+N})=Q_{n+1}(x_{k+N1},u_{k+N1})=\mathrm{}=Q_{n+N1}(x_{k+1},u_{k+1})=Q_{n+N}(x_k,u_k)`$

where $`x_{k\pm N}=x_k,u_{k\pm N}=u_k`$. Therefore the solutions turn out to be

$$\{\begin{array}{cc}\{\begin{array}{c}\overline{x}_k=x_{k1}\frac{Q_N(x_k,u_k)}{Q_1(x_k,u_k)}\hfill \\ \overline{u}_k=u_k\frac{Q_1(x_k,u_k)}{Q_N(x_k,u_k)}\hfill \end{array}\hfill & \text{A-type},\\ \{\begin{array}{c}\overline{x}_k=u_k\hfill \\ \overline{u}_k=x_{k1}\hfill \end{array}\hfill & \text{B-type}.\end{array}$$
(83)
Similarly the time reversal map is obtained by writing eq.(77) as

$`\{\begin{array}{ccc}\begin{array}{ccc}x_k\hfill & =& {\displaystyle \frac{(\overline{x}_k+\overline{u}_{k+1})x_{k1}\overline{x}_k\overline{u}_k}{x_{k1}}}\hfill \\ u_k\hfill & =& {\displaystyle \frac{(\overline{x}_k+\overline{u}_{k+1})u_{k+1}\overline{x}_{k+1}\overline{u}_{k+1}}{u_{k+1}}}\hfill \end{array}\hfill & & \end{array}`$ (87)
corresponding to the map $`(x_{k1},u_{k+1})(x_k,u_k)`$.

This is again a modified 5th pattern of the Möbius transformation. The modification matrix and the elements of the map are given by

$$(\beta _n,\gamma _n)=\{\begin{array}{c}(\overline{x}_{k+1n},\overline{u}_{k+2n})\hfill \\ (\overline{u}_{k+n},\overline{x}_{k1+n})\hfill \end{array},U_n=\left(\begin{array}{cc}1& 0\\ 0& \gamma _n\end{array}\right),$$
(88)
and the solutions

$$\{\begin{array}{cc}\{\begin{array}{c}x_{k1}=\overline{x}_k\frac{Q_N^{}(\overline{x}_k,\overline{u}_{k+1})}{Q_1^{}(\overline{x}_k,\overline{u}_{k+1})}\hfill \\ u_k=\overline{u}_k\frac{Q_1^{}(\overline{x}_k,\overline{u}_{k+1})}{Q_N^{}(\overline{x}_k,\overline{u}_{k+1})}\hfill \end{array}\hfill & \text{A-type},\\ \{\begin{array}{c}x_k=\overline{u}_{k+1}\hfill \\ u_k=\overline{x}_k\hfill \end{array}\hfill & \text{B-type}.\end{array}$$
(89)
Instead of eq.(80) and eq.(87) we could write the dToda equation eq.(77) in the form

$`\{\begin{array}{ccc}\overline{x}_k\hfill & =& {\displaystyle \frac{x_{k1}u_k}{\overline{x}_{k1}(x_{k1}+u_{k1})}}\hfill \\ \overline{u}_k\hfill & =& {\displaystyle \frac{x_{k1}u_k}{\overline{u}_{k+1}(x_k+u_k)}}\hfill \end{array},`$ (92)
$`\{\begin{array}{ccc}x_k\hfill & =& {\displaystyle \frac{\overline{x}_{k+1}\overline{u}_{k+1}}{x_{k+1}(\overline{x}_{k+1}+\overline{u}_{k+1})}}\hfill \\ u_k\hfill & =& {\displaystyle \frac{\overline{x}_k\overline{u}_k}{u_{k1}(\overline{x}_{k1}+\overline{u}_k)}}\hfill \end{array}`$ (95)
corresponding to $`(\overline{x}_{k1},\overline{u}_{k+1})(\overline{x}_k,\overline{u}_k)`$ and $`(x_{k+1},u_{k1})(x_k,u_k)`$. These are modified 6th patterns whose elements and the modification matrix are given by

$$(\beta _n,\gamma _n)=\{\begin{array}{cc}\{\begin{array}{c}(x_{k1n},u_{k1n})\hfill \\ (u_{k+n},x_{k+n})\hfill \end{array}\hfill & \\ \{\begin{array}{c}(\overline{x}_{k+1+n},\overline{u}_{k+1+n})\hfill \\ (\overline{u}_{kn},\overline{x}_{k1n})\hfill \end{array}\hfill & \end{array},U_n=\left(\begin{array}{cc}1/\gamma _{n1}& 0\\ 0& 1\end{array}\right)$$
(96)
for eq.(92) and eq.(95), respectively. Applying these data to eq.(68) we obtain the same results of eq.(83) and eq.(89).

### 3.2 Generalizations

#### 3.2.1 Symmetry of B-type solutions

We have seen that the B-type solution is simply 1 if the map $`M_n`$ is the form of either 5th or 6th pattern. A modification of the matrix by $`U_n`$ will change the solutions as well as the equation while the discriminant remains constant. In particular a different B-type solution will correspond to a different set of equations, and hence different A-type solution. If the B-type solution is rational the A-type solution must be also rational satisfying the same equations.

To see how it works we try to generalize the B-type solution eq.(83) to

$$\{\begin{array}{c}\overline{x}_k=u_{k+i}\hfill \\ \overline{u}_k=x_{k+j}\hfill \end{array},i,j=0,1,2,\mathrm{},N1.$$
(97)
It is not difficult to find their corresponding equations

$$\{\begin{array}{c}\overline{x}_{k+l}\overline{u}_k=x_{k+j}u_{k+i+l}\hfill \\ \overline{x}_k+\overline{u}_{k+m}=x_{k+j+m}+u_{k+i}\hfill \end{array},l,m=0,1,2,\mathrm{},N1.$$
(98)
From the construction this equation has also A-type solution which is a rational polynomial of $`x`$ and $`u`$. And we can set $`h=m+l=\pm 1`$ without loss of generality.

We would like to emphasize that the B-type solution follows automatically to the symmetry possessed by the equation. In other words, the rational solution of A-type is a consequence of the symmetry of the Toda lattice equation.

#### 3.2.2 $`q`$-difference version

In the paper  Kajiwara et al. have studied affine-Weyl symmetry of the discrete time Painlevé IV equation with parameters, which is called $`q`$P<sub>IV</sub>. We have shown in our previous paper  that $`q`$P<sub>IV</sub> is one of the equations which are characterized by the square root free maps.

Now we want to know whether exist similar generalization of the discrete time Toda lattice eq.(77) which we have been studying. To this end we consider the following equation with a set of parameters $`a_k,b_k,k=1,2,\mathrm{},N`$ :

$$\{\begin{array}{ccc}\frac{\overline{x}_{k+l}}{a_{k+l}}\frac{\overline{u}_k}{b_k}& =& a_{k+j}x_{k+j}b_{k+i+l}u_{k+i+l}\\ \frac{\overline{x}_k}{a_k}+\frac{\overline{u}_{k+m}}{b_{k+m}}& =& a_{k+j+m}x_{k+j+m}+b_{k+i}u_{k+i}.\end{array}$$
(99)
Let us call this equation a $`q`$-difference generalized discrete time Toda equation ($`q`$Toda). Remarkably this equation has the symmetry of modified 5th or 6th pattern of the Möbius map. So, this is transformed to a quadratic equation of every $`\overline{x},\overline{u},x`$ or $`u`$, and has two types of map as solutions to eq.(99).

In fact if we introduce

$`M_n`$ $`=`$ $`\left(\begin{array}{cc}a_{k+s+nh}x_{k+s+nh}+b_{k+i+nh}u_{k+i+nh}& a_{k+s+nh}x_{k+s+nh}\\ b_{k+i+nh}u_{k+i+nh}& 0\end{array}\right),`$ (102)
$`U_n`$ $`=`$ $`\left(\begin{array}{cc}1& 0\\ 0& b_{k+i+nh}u_{k+i+nh}\end{array}\right),`$ (105)
we obtain two maps as follows :

$$\frac{\overline{x}_k}{a_k}=\{\begin{array}{cc}a_{k+sh}x_{k+sh}\frac{Q_N^h(a_{k+s}x_{k+s},b_{k+i}u_{k+i})}{Q_1^h(a_{k+s}x_{k+s},b_{k+i}u_{k+i})},\hfill & \text{A-type}\hfill \\ & \\ b_{k+i}u_{k+i},\hfill & \text{B-type}\hfill \end{array}$$
(106)
where $`s=j+m,h=l+m`$, and $`Q_n^h`$ is defined by

$$Q_n^h(x_k,u_k):=\{\begin{array}{ccc}Q_n(x_k,u_k)\hfill & \text{if}& h=1\hfill \\ Q_n^{}(x_k,u_k)\hfill & \text{if}& h=1.\hfill \end{array}$$

#### 3.2.3 Higher dimension

We can further generalize $`q`$Toda to higher dimensions by increasing the number of suffices:

$$(x_k,u_k)(x_{k_1,k_2,\mathrm{},k_M},u_{k_1,k_2,\mathrm{},k_M}).$$

The B-type map should have the form

$$\{\begin{array}{ccc}\frac{\overline{x}_{k_1,k_2,\mathrm{},k_M}}{a_{k_1,k_2,\mathrm{},k_M}}& =& b_{k_1+i_1,k_2+i_2,\mathrm{},k_M+i_M}u_{k_1+i_1,k_2+i_2,\mathrm{},k_M+i_M}\\ & & \\ \frac{\overline{u}_{k_1,k_2,\mathrm{},k_M}}{b_{k_1,k_2,\mathrm{},k_M}}& =& a_{k_1+j_1,k_2+j_2,\mathrm{},k_M+j_M}x_{k_1+j_1,k_2+j_2,\mathrm{},k_M+j_M}.\end{array}$$
(107)
In the case of $`M=2`$, for instance,

$$\{\begin{array}{ccc}\frac{\overline{x}_{k,l}}{a_{k,l}}& =& b_{k+1,l+1}u_{k+1,l+1}\\ & & \\ \frac{\overline{u}_{k,l}}{b_{k,l}}& =& a_{k+1,l1}x_{k+1,l1}\end{array}$$
(108)
satisfies

$$\{\begin{array}{ccc}\frac{\overline{x}_{k1,l}\overline{u}_{k1,l}}{a_{k1,l}b_{k1,l}}& =& a_{k,l1}b_{k,l+1}x_{k,l1}u_{k,l+1},\\ & & \\ \frac{\overline{x}_{k,l}}{a_{k,l}}+\frac{\overline{u}_{k1,l+1}}{b_{k1,l+1}}& =& a_{k,l}x_{k,l}+b_{k+1,l+1}u_{k+1,l+1}.\end{array}$$
(109)
This is the $`q`$KP equation which we derived in our previous paper.

### 3.3 Reductions

It has been well known that the discrete time Lotka-Volterra equation (dLV) can be obtained from the discrete time Toda equation by a Miura transformation . We like to show, in this section, that the same dLV-type as well as the discrete time KdV (dKdV)-type equations can be derived from the discrete time Toda lattice eq.(99) by simple reductions. The basic idea of our reduction is to deal with one of two equations in $`q`$Toda as a constraint to eliminate half of the variables.

#### 3.3.1 Reduction to the dLV-type

The way of reduction to a dLV-type equation is to impose an additional constraint to the second equation of $`q`$Toda as follows.

$$\{\begin{array}{ccccc}\frac{\overline{x}_{k+l}}{a_{k+l}}\frac{\overline{u}_k}{b_k}& =& a_{k+j}x_{k+j}b_{k+i+l}u_{k+i+l}& & \\ \frac{\overline{x}_k}{a_k}+\frac{\overline{u}_{k+m}}{b_{k+m}}& =& a_{k+j+m}x_{k+j+m}+b_{k+i}u_{k+i}& =& r_k\end{array}$$
(110)
where $`r_k`$’s are constants.

This constraint enables us to eliminate variables, say $`u`$ and $`\overline{u}`$, from the first equation. Namely the following substitution

$$\{\begin{array}{ccc}\frac{\overline{u}_{k+m}}{b_{k+m}}& & \frac{\overline{x}_k}{a_k}r_k\\ b_{k+i}u_{k+i}& & a_{k+j+m}x_{k+j+m}r_k\end{array}$$
(111)
yields the dLV-type equation

$$\frac{\overline{x}_{k+l}}{a_{k+l}}\left(r_{km}+\frac{\overline{x}_{km}}{a_{km}}\right)=a_{k+j}x_{k+j}\left(r_{k+l}+a_{k+j+m+l}x_{k+j+m+l}\right).$$
(112)
When $`r_k=a_k=b_k=1,k=1,2,\mathrm{},N`$, this is the dLV equation. Various types of $`q`$LV equations which were discussed in our previous paper, are included in this equation.

It is straightforward to obtain two maps of this equation. We find

$$\frac{\overline{x}_k}{a_k}=\{\begin{array}{cc}a_{k+sh}x_{k+sh}\frac{Q_N^h(a_{k+s}x_{k+s},r_ka_{k+j+m}x_{k+j+m})}{Q_1^h(a_{k+s}x_{k+s},r_ka_{k+j+m}x_{k+j+m})},\hfill & \text{A-type}\hfill \\ & \\ r_ka_{k+j+m}x_{k+j+m},\hfill & \text{B-type}\hfill \end{array}$$
(113)
which could be also obtained from the map of qToda eq.(99) by the substitution eq.(111). The $`q`$P<sub>IV</sub> of Kajiwara, Noumi, Yamada  is included as one of the A-type map.

#### 3.3.2 Reduction to the dKdV-type

Let us try an alternative case of constraint. Namely we impose a constraint to the first equation of eq.(99).

$$\{\begin{array}{ccccc}\frac{\overline{x}_{k+l}}{a_{k+l}}\frac{\overline{u}_k}{b_k}& =& a_{k+j}x_{k+j}b_{k+i+l}u_{k+i+l}& =& r_k,\\ \frac{\overline{x}_k}{a_k}+\frac{\overline{u}_{k+m}}{b_{k+m}}& =& a_{k+j+m}x_{k+j+m}+b_{k+i}u_{k+i}.& & \end{array}$$
(114)
This enables us to eliminate $`u`$ and $`\overline{u}`$ from the second equation according to the rule

$$\{\begin{array}{ccc}\frac{\overline{u}_k}{b_k}& & \frac{a_{k+l}r_k}{\overline{x}_{k+l}}\\ b_{k+i+l}u_{k+i+l}& & \frac{r_k}{a_{k+j}x_{k+j}}\end{array}$$
(115)
from which we obtain

$$\frac{\overline{x}_k}{a_k}+\frac{a_{k+l+m}r_{k+m}}{\overline{x}_{k+l+m}}=a_{k+j+m}x_{k+j+m}+\frac{r_{kl}}{a_{k+jl}x_{k+jl}}.$$
(116)
This equation containes the discrete time KdV equation as a special case. Hence we call eq.(116) the $`q`$KdV equation. The corresponding solutions are given again by simply replacing $`u`$ and $`\overline{u}`$ in the solutions of $`q`$Toda according to eq.(115).

$$\frac{\overline{x}_k}{a_k}=\{\begin{array}{cc}a_{k+sh}x_{k+sh}\frac{Q_N^h(a_{k+s}x_{k+s},r_{kl}/a_{k+jl}x_{k+jl})}{Q_1^h(a_{k+s}x_{k+s},r_{kl}/a_{k+jl}x_{k+jl})},\hfill & \text{A-type}\hfill \\ & \\ r_{kl}/a_{k+jl}x_{k+jl}.\hfill & \text{B-type}\hfill \end{array}$$
(117)
### 3.4 Reflection maps

We define some symbols of the two types of map and their inverse maps of dToda as follows,

$$\begin{array}{ccccccccccc}𝒜\hfill & :& \left(\begin{array}{c}x_k\\ u_k\end{array}\right)\hfill & & \left(\begin{array}{c}x_{k1}\frac{Q_N(x_k,u_k)}{Q_1(x_k,u_k)}\\ u_k\frac{Q_1(x_k,u_k)}{Q_N(x_k,u_k)}\end{array}\right),\hfill & & \hfill & :& \left(\begin{array}{c}x_k\\ u_k\end{array}\right)\hfill & & \left(\begin{array}{c}u_k\\ x_{k1}\end{array}\right),\hfill \\ 𝒜^1\hfill & :& \left(\begin{array}{c}x_{k1}\\ u_k\end{array}\right)\hfill & & \left(\begin{array}{c}x_k\frac{Q_N^{}(x_k,u_{k+1})}{Q_1^{}(x_k,u_{k+1})}\hfill \\ u_k\frac{Q_1^{}(x_k,u_{k+1})}{Q_N^{}(x_k,u_{k+1})}\hfill \end{array}\right),\hfill & & ^1\hfill & :& \left(\begin{array}{c}x_k\\ u_k\end{array}\right)\hfill & & \left(\begin{array}{c}u_{k+1}\hfill \\ x_k\hfill \end{array}\right).\hfill \end{array}$$

These maps satisfy the following relations.

$$\begin{array}{cc}(\mathrm{I})& 𝒜𝒜^1=𝒜^1𝒜=1,^1=^1=1,\hfill \\ (\mathrm{II})& 𝒜^2=^2𝒜,\hfill \\ (\mathrm{III})& 𝒜𝒜=^3,\hfill \end{array}$$
(118)
or

$$\begin{array}{cc}(\mathrm{I}^{})& ^1=^1=1,\hfill \\ (\mathrm{II}^{})& 𝒞=𝒟,\hfill \\ (\mathrm{III}^{})& 𝒞^2=𝒟^2=1,\hfill \end{array}$$
(119)
where $`𝒞=𝒜^1=^1𝒜`$ and $`𝒟=𝒜^1=𝒜^1`$. The relations eq.(119) are derived from eq.(118) by some calculations, automatically.

The reflection maps $`𝒞,𝒟`$ change the inverse A-type map to the inverse B-type one and the A-type one to the B-type one, respectively. Namely,

$$\begin{array}{cccc}𝒞:& \text{inverse A-type}\hfill & & \text{inverse B-type},\hfill \\ & \text{inverse B-type}\hfill & & \text{inverse A-type},\hfill \\ 𝒟:& \text{A-type}\hfill & & \text{B-type},\hfill \\ & \text{B-type}\hfill & & \text{A-type}.\hfill \end{array}$$

Hence, these reflection maps are permutations of solutions of the quadratic equation eq.(14) in the case of dToda, i.e. Galois tranformations of eq.(14). Furthermore, the generalized version of dToda eq.(98) satisfy these relations eq.(118) or eq.(119), but $`q`$Toda does not. These reflection maps play an important role in affine-Wyel symmetries and $`q`$-difference Painlevé equations as studied by Kajiwara, et al. .

## 4 Discrete Time Evolution of the $`q`$Toda

A sequence of B-type map does not generate an orbit but causes exchange of the variables $`x`$ and $`u`$, whereas an A-type map generates an orbit. We will study a discrete time evolution of A-type map in this section. To make the time dependence explicit we write the variables $`\overline{x}`$ and $`\overline{u}`$ as $`x^{t+1}`$ and $`u^{t+1}`$ in the following.

### 4.1 Time evolution of the A-type maps

A repeated use of eq.(83) enables one to express $`x_k^{t+1}`$ and $`u_k^{t+1}`$ in terms of $`x_k^0`$’s and $`u_k^0`$’s as

$$\{\begin{array}{c}x_k^{t+1}=x_{k1}^t\frac{Q_{k1}^t}{Q_k^t}=\mathrm{}=x_{kt1}^0\underset{s=0}{\overset{t}{}}\frac{Q_{ks1}^{ts}}{Q_{ks}^{ts}},\hfill \\ u_k^{t+1}=u_k^t\frac{Q_k^t}{Q_{k1}^t}=\mathrm{}=u_k^0\underset{s=0}{\overset{t}{}}\frac{Q_k^{ts}}{Q_{k1}^{ts}}.\hfill \end{array}$$
(120)
Here we denote by $`Q_k^t`$ the following polynomial of $`x^t`$’s and $`u^t`$’s,

$$Q_k^t:=\underset{j=1}{\overset{N}{}}u_{k+1}^tu_{k+2}^t\mathrm{}u_{k+j1}^tx_{k+j}^t\mathrm{}x_{k+N1}^t.$$
(121)
If we substitute eq.(120) into the right hand side of eq.(121) we will find a considerable cancellation of factors in the numerator and the denominator. The cancellation might be associated with a reduction of algebraic entropy . As a consequence we obtain a recurrence equation for $`Q^t`$:

$$Q_k^{t+1}=\underset{j=1}{\overset{N}{}}u_{k+1}^0u_{k+2}^0\mathrm{}u_{k+j1}^0x_{k+jt1}^0\mathrm{}x_{k+Nt2}^0\underset{s=0}{\overset{t}{}}\frac{Q_{k+j1}^{ts}Q_{k+js1}^{ts}}{Q_k^{ts}Q_{k+Ns1}^{ts}}.$$
(122)
The A-type solutions are given if $`Q^t`$ is found by solving eq.(122) iteratively. The generalization of this result to the case of $`q`$Toda is straightforward. This is, however, a cumbersome task. And so, we will look orbits of A-type map described by numerical computation.

### 4.2 Numerical observation of ultra-discrete version of the A-type maps

It is not difficult to see behaviour of orbits of the A-type map eq.(120) as long as we generate them numerically by using computers. They are constrained on curves since there are sufficient number of conserved quantities. If parameters $`a_k`$’s and $`b_k`$’s are introduced the orbits tend to expand and closed curves turn to spirals. These are behaviours which we could expect naively from our previous experiences. Instead of studying them in detail we would like to present here behaviour of ultra-discrete version of solutions.

It has been known that if a discrete integrable equation has its ultra-discrete version the integrability is preserved . There exists a systematic method to ultra-discretize a discrete integrable equation if it contains only product and summation of variables but not subtractions. Namely we introduce a new variable $`X`$ and replace $`x`$ by $`x=e^{X/ϵ}`$. The ultra-discretization procedure  is defined by

$$\underset{ϵ+0}{lim}\left(e^{X/ϵ}e^{Y/ϵ}\right)X+Y,\underset{ϵ+0}{lim}\left(e^{X/ϵ}+e^{Y/ϵ}\right)\mathrm{max}(X,Y).$$

To be specific let us consider the case of dKdV-type equation. The $`q`$-difference version has been given by eq.(116). If we further specify to 3-point case the equations are

$$\frac{\overline{x}_k}{a_k}+\frac{a_{k+1}}{\overline{x}_{k+1}}=a_{k+1+j}x_{k+1+j}+\frac{1}{a_{k+j}x_{k+j}},$$
(123)
$$k=1,2,3,j=0,1,2,a_{k+3}=a_k,x_{k+3}=x_k.$$

The corresponding A-type map can be also found readily from the one of $`q`$-Toda equation by the reduction eq.(115).

$$\overline{x}_{kj}=a_{kj}a_{k+2}x_{k+2}\frac{1+a_{k+2}a_{k+1}x_{k+2}x_{k+1}+a_ka_{k+1}^2a_{k+2}x_kx_{k+1}^2x_{k+2}}{1+a_ka_{k+2}x_kx_{k+2}+a_ka_{k+1}a_{k+2}^2x_kx_{k+1}x_{k+2}^2},$$
(124)
$$k=1,2,3,j=0,1,2,a_{k+3}=a_k,x_{k+3}=x_k.$$

The ultra-discrete version of eq.(124) turns to

$$\begin{array}{cc}\overline{X}_{kj}=& A_{kj}+A_{k+2}+X_{k+2}\hfill \\ & +\mathrm{max}(0,A_{k+2}+A_{k+1}+X_{k+2}+X_{k+1},A_k+2A_{k+1}+A_{k+2}+X_k+2X_{k+1}+X_{k+2})\hfill \\ & \mathrm{max}(0,A_k+A_{k+2}+X_k+X_{k+2},A_k+A_{k+1}+2A_{k+2}+X_k+X_{k+1}+2X_{k+2}).\hfill \end{array}$$
(125)
Here we set $`a_k=e^{A_k/ϵ},k=1,2,3`$. Let $`X_k`$’s and $`A_k`$’s be rational numbers and impose a condition $`A_1+A_2+A_3=0`$ so that the map returns to initial values after some steps. Under these circumstances we find maps such as presented in Figures 2 to 4.

Figure 1 shows the case of $`A_k=0,k=1,2,3`$ for all j=0,1,2. Similarly Fig. 2 corresponds to the case of $`A_1+A_2+A_3=0`$ and $`j=0`$. The straight lines in these figures are known to correspond to integrable orbits . Figures 3 and 4 are the cases of $`j=1`$ and $`j=2`$ respectively. They do not seem to describe integrable orbits. Therefore the orbits behave differently for different value of $`j`$. This must be contrasted with the case of dKdV in which orbits form smooth closed curves under the condition $`a_1a_2a_3=1`$ for all $`j=0,1,2`$.

It will be worth while to point out that the dKdV map eq.(124) in the case of $`j=0`$ is invariant under the following Bäcklund transformation

$$\begin{array}{cc}\pi (a_k)=a_{k+1},\hfill & \pi (x_k)=x_{k+1},\hfill \\ s_i(a_j)=a_j(a_{i+1}a_{i1})^{\lambda _{ij}},\hfill & s_i(x_j)=x_j\left(\frac{a_{i+1}a_{i1}+x_{i+1}x_{i1}}{1+a_{i+1}a_{i1}x_{i+1}x_{i1}}\right)^{\gamma _{ij}},\hfill \end{array}$$

where

$$\mathrm{\Lambda }:=(\lambda _{ij})=\left(\begin{array}{ccc}\hfill 2& \hfill 1& \hfill 1\\ \hfill 1& \hfill 2& \hfill 1\\ \hfill 1& \hfill 1& \hfill 2\end{array}\right),\mathrm{\Gamma }:=(\gamma _{ij})=\left(\begin{array}{ccc}\hfill 0& \hfill 1& \hfill 1\\ \hfill 1& \hfill 0& \hfill 1\\ \hfill 1& \hfill 1& \hfill 0\end{array}\right),$$

$$s_k^2=1,\pi ^3=1,s_ks_{k+1}s_k=s_{k+1}s_ks_{k+1},\pi s_k=s_{k+1}\pi $$

and $`i,j,k=1,2,3`$ and $`a_{k+3}=a_k,x_{k+3}=x_k`$. This is the same affine Weyl symmetry shared by $`q`$P<sub>IV</sub> . But $`N`$($``$4)-point dKdV maps of $`j=0`$ does not have this symmetry even though it exhibits an integrable orbit.

We have also studied 3-point ultra-discrete LV-type maps in the same manner. The results are similar to those of KdV-type maps and found one integrable map corresponding to the $`q`$P<sub>IV</sub>.

## 5 Summary

The main purpose of this paper was to clarify the method characterizing discrete integrable systems, which we introduced in our porevious works. We have shown that, in the case of integrable systems, two adjacent variables are related by a Möbius transformation of special forms. We found at least 6 patterns of such forms.

This observation enabled us to obtain quite large class of new discrete integrable systems. The generalized discrete Toda equations (dToda) of eq.(98) is one of such examples. By certain reductions of the dToda equation dKdV and dLV were obtained. The method admits introduction of arbitrary parameters into the equations and further generalization to $`q`$Toda systems.

If we solve the equations in the form in which each variable is given as a function of variables in the previous time, there exist always two types of solutions which we called A-type and B-type. The B-type map does not create orbit but exchanges variables, while a sequence of the A-type map generates an orbit.

The map of A-type is written as a rational polynomial of variables whose numerator and denominator are expressed by minor determinants of the matrix $`Q`$ in eq.(69). It is interesting to notice that this matrix $`Q`$ is the same one known in the field of affine Weyl group and $`q`$-Painlevé equations .

We have not succeeded analyzing full behaviour of a sequence of A-type maps. It is one of the most interesting problems to find behaviour of the map when parameters $`a_k`$’s and $`b_k`$’s are introduced. Some numerical results using ultra-discrete method have shown that the introduction of the parameters preserve integrability in some cases but not always. This result seems to show that the non square root map condition is not sufficient but necessary for a discrete map being integrable. It is desirable to clarify this point further.


format: text

text:
# Line formation in convective stellar atmospheres

## 1 Introduction

Theoretical model atmospheres serve as the fundamental tool for the quantitative analysis of stellar spectra. The basic parameters of a stellar atmosphere – effective temperature, surface gravity, and chemical composition – can be inferred from the comparison between observed and calculated spectra. Much effort has been spent to improve the theoretical models by incorporating as much realistic physics as possible. Specifically, in late-type stars one is confronted with the problem of a detailed treatment of both radiative and convective energy transport. A high degree of sophistication has been achieved in modeling radiative transfer: standard stellar atmosphere codes can now handle the influence of millions of spectral lines.

In contrast, the treatment of convection is still rather crude. Usually, stellar atmospheres are described as 1-dimensional hydrostatic configurations where the convective energy transport is calculated from the so-called mixing-length theory (MLT, Vitense V53 (1953); Böhm-Vitense BV58 (1958)) or variants thereof (e.g. Canuto et al.CGM96 (1996)). However, modern spectroscopic observations have reached a level of precision that requires improved model atmospheres for a reliable interpretation, including a proper treatment of hydrodynamical phenomena.

From the point of view of standard model atmospheres, convection affects the *temperature structure* in a twofold way: it influences the *mean vertical stratification* and introduces *horizontal inhomogeneities*. MLT is designed to model only the mean structure and is not suitable to construct realistic multi-component model atmospheres. Although realistic hydrodynamical model atmospheres do exist for the Sun and a few solar-type stars due to the pioneering work by Nordlund (No82 (1982)), Nordlund & Dravins (ND90 (1990)), Stein & Nordlund (SN89 (1989), SN98 (1998)), and Ludwig et al. (LFS99 (1999)), to name just a few important contributions, such complex models are not usually applied for the analysis of stellar spectra.

Observationally, the solar granulation is the visible imprint of thermal convection extending far into the photospheric layers where spectral lines are formed. The effect of the associated photospheric temperature inhomogeneities on the line formation process and the consequences for spectroscopic abundance determinations are not yet well understood, but it is generally believed that the errors introduced by representing a dynamic, inhomogeneous atmosphere by a static, plane-parallel model are small. However, reliable quantitative estimates of this effect are difficult to find. Earlier studies suffer from uncertainties in the underlying empirical two-component models, which are rather *ad hoc* and lack a firm physical foundation (e.g. Hermsen He82 (1982)). Later investigations based on relatively crude numerical convection models (Atroshchenko & Gadun AG94 (1994); Gadun & Pavlenko GP97 (1997)) led to questionable and partly ambiguous conclusions. Recent determinations of the solar photospheric Fe and Si abundance by Asplund et al. (2000b ) and Asplund (2000a ), respectively, rely on state-of-the-art 3D numerical convection simulations and take all hydrodynamical effects fully into account. However, comparison of these results with those based on the empirical 1D Holweger-Müller (HM74 (1974)) solar atmosphere does not allow to separate the influence of inhomogeneities from the effect of the mean stratification, and to study the 1D/3D abundance corrections as a function of the line parameters in a systematic way.

In this work we investigate the *spectroscopic effects of horizontal temperature inhomogeneities in the solar atmosphere* based on detailed 2D hydrodynamical models of solar surface convection (Freytag et al. FLS96 (1996); Ludwig et al. LFS99 (1999), Steffen S00 (2000)). We do *not* address the question whether the mean temperature structure supplied by MLT is an appropriate representation of the mean stratification obtained from hydrodynamical simulations (the problem of the ‘right’ choice of the mixing-length parameter $`\alpha `$ for recovering the correct mean temperature stratification was investigated by Steffen et al. (SLF95 (1995)) in the context of white dwarf atmospheres). Rather, our main question is: how large are the systematic errors of standard spectroscopic abundance determinations due to the fact that one replaces the ‘real’ line profile, which actually is the result of averaging the spatially resolved line *intensities* over the granulation pattern, by the line profile computed from a single plane-parallel temperature stratification?

In Sect. 2, we introduce a simple analytical argument to illustrate the basic effect of line strengthening in the presence of temperature fluctuations. Then we summarize the main characteristics of the numerical simulations of solar surface convection in Sect. 3, outline our method to derive ‘granulation abundance corrections’ from a differential 1D/2D comparison of computed line equivalent widths in Sect. 4, before finally presenting the results for the Sun in Sect. 5. A summary of our findings is given in Sect. 6.

In order to give a physical explanation for the numerical results obtained in this work, we have carried out a detailed investigation the process of LTE line formation in an inhomogeneous stellar atmosphere based on the analysis of the transfer equation and the evaluation of line depression contribution functions. This study is presented in the second paper of this series (Paper II).

## 2 A simple analytical illustration

The strength of a spectral line formed in a plane-parallel (1D) atmosphere will differ from the strength of the same line formed in an inhomogeneous atmosphere with the same chemical composition and mean temperature stratification because, in general, there is a *non-linear relationship between line strength and temperature*. It may be expected that lines with a very high excitation potential are most affected because (i) they originate from deep layers where convection is well developed, and (ii) they are highly sensitive to temperature fluctuations. Likewise, the molecular equilibrium is known to be highly sensitive to temperature fluctuations, so molecular lines might also be strongly affected if the temperature inhomogeneities reach the higher photosphere where these lines are formed.

As an illustration of the expected effect, consider the following simple analytical argument. Assume that the equivalent width of a weak line $`W`$ depends exponentially on the inverse temperature $`\theta =5040/T`$, as suggested by the form of the Saha-Boltzmann equations, 111Note that Eq.(1) is a reasonable approximation only for weak lines of ‘majority species’ with high excitation potential (cf. Gray Gr92 (1992)).

$$W(\mathrm{\Delta }\theta )=W_0\mathrm{exp}\left\{E_i\mathrm{\Delta }\theta \mathrm{ln}10\right\},$$
(1)
where $`E_i`$ is the excitation potential in eV and $`\mathrm{\Delta }\theta =\theta \theta _0`$. Assume further that $`\theta `$ is normally distributed about the mean $`\theta _0`$ with standard deviation $`\sigma `$,

$$g(\mathrm{\Delta }\theta )=\frac{1}{\sigma \sqrt{2\pi }}\mathrm{exp}\left\{\frac{\mathrm{\Delta }\theta ^2}{2\sigma ^2}\right\}.$$
(2)
Then it can be shown by analytical integration that

$`W`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}g(\mathrm{\Delta }\theta )W(\mathrm{\Delta }\theta )d\mathrm{\Delta }\theta `$ (3)
$`=`$ $`W_0\mathrm{exp}\left\{{\displaystyle \frac{1}{2}}E_i^2\sigma ^2(\mathrm{ln}10)^2\right\}`$
i.e. the mean value $`W`$ of the equivalent width formed in an inhomogeneous medium is always greater than the equivalent width for the unperturbed medium, $`W_0`$. As long as $`E_i\sigma `$ is small enough, the effect is negligible, but this is not true in general. Assuming $`\sigma =0.05`$, corresponding to 6% temperature fluctuation around $`T_0`$ = 6000 K, and $`E_i=10`$ eV, the effect is very substantial: $`W2W_0`$, i.e. the elemental abundance derived from such lines would be *overestimated* by a factor of 2! Note that, for example, several important neutral oxygen lines have $`E_i`$ between 9 and 11 eV.

A similar result is obtained when assuming a totally different distribution of $`\mathrm{\Delta }\theta `$,

$$f(\mathrm{\Delta }\theta )=\frac{1}{\pi \sqrt{2\sigma ^2\mathrm{\Delta }\theta ^2}},$$
(4)
corresponding to sinusoidal temperature fluctuations in one horizontal direction,

$$\mathrm{\Delta }\theta =\sigma \sqrt{2}\mathrm{sin}(2\pi x),$$
(5)
where $`x`$ is the spatial coordinate. According to Eq.(1), the fluctuation of the equivalent width is then

$$W(x)=W_0\mathrm{exp}\left\{E_i\mathrm{ln}10\sigma \sqrt{2}\mathrm{sin}(2\pi x)\right\}.$$
(6)
Adopting again $`\sigma =0.05`$ and $`E_i=10`$ eV, we have plotted $`W(x)`$ in Fig. 1. Clearly, the fluctuation of $`W`$ is highly asymmetric with respect to $`W_0`$, and numerical averaging over one spatial wavelength gives $`W1.8W_0`$.

Looking at Table 1, it turns out that the simple model severely overestimates the abundance correction in the case of the oxygen lines mentioned as an example above. This is because, in reality, the line formation process is much more complicated, since the line strength depends not only on the line opacity but also on the continuum opacity and on the local temperature gradient of the atmosphere (see Paper II). Moreover, spectral lines usually form over a substantial depth range, and the sign of $`\mathrm{\Delta }\theta `$ may depend on depth, leading to partial cancellation of the local effects, especially for stronger lines. Finally, saturation can play a significant role as well. Hence, realistic abundance corrections can only be derived from detailed line formation calculations based on multi-dimensional hydrodynamical model atmospheres.

## 3 Time-dependent Radiation Hydrodynamics

In order to adequately address the complex problem outlined above, we have investigated the effects of temperature inhomogeneities on the formation of a variety of spectral lines by means of a realistic 2-D radiation hydrodynamics (RHD) simulation of solar surface convection. In the following, we briefly summarize the foundations of the hydrodynamical simulation and the main results concerning the solar photosphere.

### 3.1 2-dimensional numerical simulation

Stellar surface convection is governed by the conservation equations of hydrodynamics, coupled with the equations of radiative energy transfer. Our hydrodynamical models result from the numerical integration of this set of partial differential equations. This approach constitutes an increasingly powerful tool to study in detail the time-dependent hydrodynamical properties of the solar granulation.

Just like ‘classical’ stellar atmospheres, the hydrodynamical models are characterized by effective temperature, $`T_{\mathrm{eff}}`$, surface gravity, $`g`$, and chemical composition of the stellar matter. But in contrast to the mixing-length models, they account for ‘overshoot’ and no longer have any free parameter to adjust the efficiency of the convective energy transport. Based on first principles, RHD simulations provide physically consistent *ab initio* models of stellar convection which can serve to address a variety of questions, including the problem of line formation in inhomogeneous stellar atmospheres.

The models used for the present investigation comprise a small section near the solar surface, extending over 7 pressure scale heights in the vertical direction, including the photosphere, the thermal boundary layer near optical depth $`\tau _{\mathrm{Ross}}=1`$, and parts of the subphotospheric layers. Only the uppermost layers of the deep solar convection zone can be included in the model, requiring an open lower boundary. The simulations are designed to resolve the solar granulation. The effect of the smaller scales, which cannot be resolved numerically, is modeled by means of a subgrid scale viscosity (so-called Large Eddy Simulation approach). Spatial scales larger than the computational box are unaccounted for. We employ a realistic equation of state, including ionization of H, HeI, HeII and H2 molecule formation. In order to avoid problematic simplifications like the diffusion or Eddington approximation, we solve the non-local radiative transfer problem along a large number ($`N_\mathrm{x}f_{\mathrm{sub}}N_\theta N_\varphi =210428\mathrm{10\hspace{0.17em}000}`$) of rays in 5 *opacity bins*, with realistic opacities accounting adequately for the influence of spectral lines (see Nordlund No82 (1982), Ludwig Lu92 (1992)). For further details see Ludwig et al. (LJS94 (1994)), Freytag et al. (FLS96 (1996)), Ludwig et al. (LFS99 (1999)), Steffen (S00 (2000)).

### 3.2 Main results

The basic features of the numerical convection model are illustrated in Fig. 2, showing a representative snapshot from a well-relaxed simulation of the solar granulation. In the convectively unstable, subphotospheric layers, fast ($`vc_s)`$ narrow downdrafts (also called plumes or jets) stand out as the most prominent feature. They are embedded in broad ascending regions (the granules) where velocities are significantly lower. We find enormous temperature differences between the cool downflows and the hot upflows ($`\mathrm{\Delta }T_{\mathrm{max}}\mathrm{5\hspace{0.17em}000}`$ K at $`z100`$ km). Near optical depth $`\tau =1`$, efficient radiative surface cooling produces a very thin thermal boundary layer over the ascending parts of the flow, where the temperature drops sharply with height ($`\mathrm{\Delta }T/\mathrm{\Delta }z80`$ K/km, see upper panel of Fig. 3) and the gas density exhibits a local inversion. Note however, that the temperature contrast is much reduced on the optical depth scale (see bottom panel of Fig. 3).

Although convectively stable according to the Schwarzschild criterion, the photosphere ($`\tau 1`$) is by no means static. Convective flows overshooting into the stable layers from below are decelerated here and deflected sideways. This can sometimes result in transonic horizontal streams which lead to the formation of shocks in the vicinity of strong downdrafts (see vertical fronts in Fig. 2 between $`\mathrm{3\hspace{0.17em}000}<x<\mathrm{3\hspace{0.17em}500}`$ km, $`z100`$ km).

Oscillations excited by the stochastic convective motions contribute to the photospheric velocity field as well. Interestingly, the oscillation periods lie in the range 150 to 500 s for the solar simulation, in close agreement with the observed 5 minute oscillations. Controlled by the balance between dynamical cooling due to adiabatic expansion and radiative heating in the spectral lines, the mean photospheric temperature stratification is slightly cooler than in radiative equilibrium, and the resulting temperature structure is not at all plane-parallel.

According to Eq.(3), the line strengthening depends on the key quantity $`\sigma =\theta _0\delta T_{\mathrm{rms}}/T_0`$, where $`\delta T_{\mathrm{rms}}=\sqrt{(TT_0)^2}`$. We have evaluated $`\delta T_{\mathrm{rms}}`$ as a function of optical depth for our 2D numerical simulation of solar convection (Fig. 4). For the problem under investigation, temperature differences have to be taken over levels of constant optical depth, rather than at constant geometrical depth. According to the simulations, typical values in the lower photosphere are in the range $`0.04<\delta T_{\mathrm{rms}}<0.06`$. The minimum near $`\mathrm{log}\tau _{\mathrm{Ross}}=0.5`$ is related to the ‘inversion’ of the temperature fluctuations at this height.

It is worth mentioning that we have recently carried out fully 3-dimensional simulations with a new code developed by B. Freytag and M. Steffen (Freytag; Steffen et al., in preparation). We find that the depth dependence of $`\delta T_{\mathrm{rms}}(\tau )`$ is qualitatively similar in 2D and 3D. However, the amplitude of the temperature fluctuations is systematically smaller at all heights in the 3D simulation. The difference is most pronounced in the higher photosphere (see Fig. 4).

### 3.3 Accuracy of 2D and 3D convection models

The straightforward application of hydrodynamics and radiative transfer simulations is certainly the most satisfactory and consistent method to treat the line formation in a convective stellar atmosphere, although the outcome is difficult to check by independent analysis. The hydrodynamical models presented here must be regarded as somewhat preliminary for several reasons. First of all, the restriction of the flow geometry to 2 dimensions must lead to deviations from the true 3D atmospheric structure. As a consequence, the convective line shift and asymmetry predicted by the 2D granulation models is quantitatively different from that obtained with 3D models (Asplund et al. 2000a ). A further problem is the relatively poor spatial resolution of the simulations (especially in the vertical direction), which is barely sufficient to resolve the thin thermal boundary layer near $`\tau =1`$. Certainly, the Reynolds number of the simulations is many orders of magnitude below realistic values. Finally, the influence of the spectral lines on the radiative energy balance can be treated only in an approximate way. Typically, a few opacity bins are used in the simulations, as compared to thousands of frequency points in standard 1D models.

The latter problems apply to state-of-the-art 3D convection simulations as well. We note that a perfect line profile fitting does not necessarily imply a perfect model atmosphere. Another important diagnostic is the center-to-limb variation of the continuum intensity. Even the best 3D solar granulation models still fail to reproduce the observed continuum properties with high accuracy (Asplund et al. ANT99 (1999)), indicating that “the temperature structure close to the continuum forming layers could be somewhat too steep”.

In view of the above mentioned limitations, we prefer to avoid a direct application of our convection simulations. Instead, we favor the somewhat less powerful but more reliable differential approach described in the following section.

## 4 Granulation abundance corrections

### 4.1 The problem

The problem we want to address is the following. Assume you have the ‘best possible’ 1D representation of the solar atmosphere, correctly reproducing the continuum intensity distribution $`I_\lambda (\lambda )`$ and its center-to-limb variation $`I_\lambda (\mu )`$. The microturbulence parameter $`\xi _{\mathrm{mic}}`$ is adjusted such as to minimize the overall dependence of the derived abundances on line strength. The empirical Holweger-Müller (HM74 (1974)) atmosphere comes close to this ideal mean temperature structure. In the presence of temperature fluctuations, however, non-linearities in the line formation process will inevitably lead to systematic errors of abundance determinations based on individual spectral lines. The present work aims at quantifying this kind of error.

We have adopted a strictly differential procedure to quantify these errors. The basic idea is to identify our 2D hydrodynamical convection simulations with the real inhomogeneous atmosphere, and the mean structure of the simulations with the ‘best possible’ 1D representation. Note that it is not important for the hydrodynamical convection models to reproduce the mean temperature structure of the real solar atmosphere with great precision, because the related errors cancel to a first approximation in our differential comparison. Similarly, the details of the hydrodynamical velocity field are of secondary importance in the present context, since such 1D/2D differences are largely eliminated by a proper choice of the microturbulence parameter. In contrast, there is no free parameter to account for the presence photospheric *temperature inhomogeneities*. This is why we focus on the abundance corrections related to temperature fluctuations $`\delta T_{\mathrm{rms}}`$ only. The most important quantity to be provided by the convection model is therefore the magnitude and height dependence of $`\delta T_{\mathrm{rms}}`$.

In the following, we outline our method to derive LTE “granulation abundance corrections”, $`\mathrm{\Delta }_{\mathrm{gran}}`$, which must be added to the logarithmic elemental abundance derived from standard 1D model atmosphere analysis in order to compensate for the effects of temperature inhomogeneities. Note that our procedure ensures that the corrections $`\mathrm{\Delta }_{\mathrm{gran}}0`$ for vanishing horizontal *temperature* and *pressure fluctuations*.

### 4.2 The procedure

First we have selected $`S=11`$ representative snapshots from our time-dependent convection simulation (see Sect. 3). These snapshots (one of them shown in Figs. 2 and 3) all have the correct effective temperature (the horizontally averaged radiative energy flux lies within $`\pm 1\%`$ of the nominal solar flux), but differ in the number and size distribution of granules and in the continuum intensity contrast ($`\delta I_{\mathrm{rms}}`$ varies between 18.0 and 26.5% at $`\lambda \mathrm{\hspace{0.17em}5000}`$ Å). For each of the 2D snapshots ($`s=\mathrm{1\hspace{0.17em}..\hspace{0.17em}11}`$) we then constructed a *corresponding 1D model* by averaging temperature $`T`$ and gas pressure $`P`$ over the corrugated surfaces of constant Rosseland optical depth:

$$T_{\mathrm{s},1\mathrm{D}}(\tau _{\mathrm{Ross}})=T_{\mathrm{s},2\mathrm{D}}^4(\tau _{\mathrm{Ross}},x)^{1/4},$$
(7)
$$P_{\mathrm{s},1\mathrm{D}}(\tau _{\mathrm{Ross}})=P_{\mathrm{s},2\mathrm{D}}(\tau _{\mathrm{Ross}},x),$$
(8)
where $`.`$ denotes averaging over the horizontal coordinate $`x`$. Note that the 1D models constructed in this way have very nearly the same radiative flux as the underlying 2D models.

Next we compute, for given elemental abundance and atomic line parameters, the equivalent width of any given spectral line from (i) the set of 2D snapshots by horizontal averaging of the spatially resolved synthetic line profile

$$W_{\mathrm{s},2\mathrm{D}}=I_{\mathrm{c},\mathrm{s}}(x)W_{\mathrm{s},2\mathrm{D}}(x)/I_{\mathrm{c},\mathrm{s}}(x),$$
(9)
where $`I_\mathrm{c}(x)`$ is the local continuum intensity, and (ii) from the corresponding 1D models, resulting in $`W_{\mathrm{s},1\mathrm{D}}`$.

For the line formation calculations, we employed the Kiel spectrum synthesis code LINFOR under the assumption of LTE. For the 1D atmospheres, we have always adopted a depth-independent microturbulence velocity of $`\xi _{\mathrm{mic}}=1`$ km/s. For snapshots from the convection simulations we have computed synthetic spectra for two different assumptions concerning the non-thermal velocity field: in case (A) we used the depth-dependent hydrodynamical velocity field $`v(x,z)`$ obtained from the simulations, in case (B) we replaced the hydrodynamical velocity field by a depth-independent microturbulence velocity $`\xi _{\mathrm{mic}}=1`$ km/s.

As the final step, we calculate the mean equivalent widths resulting from the set of 1D and 2D line profiles, $`W_{1\mathrm{D}}=\frac{1}{S}_{s=1}^SW_{\mathrm{s},1\mathrm{D}}`$ and $`W_{2\mathrm{D}}=\frac{1}{S}_{s=1}^SW_{\mathrm{s},2\mathrm{D}}`$, respectively. For truly weak lines, the (logarithmic) “granulation abundance correction” would simply be

$$\mathrm{\Delta }_{\mathrm{gran}}^{\mathrm{weak}\mathrm{lines}}=\mathrm{log}W_{1\mathrm{D}}\mathrm{log}W2\mathrm{D}.$$
(10)
In general, however, the equivalent width is not proportional to the elemental abundance due to saturation. We account for these curve-of-growth effects in the following way. Knowing that the ‘true’ abundance, $`N_{2\mathrm{D}}`$ produces the equivalent width $`W_{2\mathrm{D}}`$ with the 2D inhomogeneous model atmosphere (and $`W_{1\mathrm{D}}`$ with the 1D mean structure), we iteratively determine the elemental abundance $`N_{1\mathrm{D}}`$ which produces the equivalent width $`W_{2\mathrm{D}}`$ with the 1D mean stratification. The “granulation abundance correction” is then given by

$$\mathrm{\Delta }_{\mathrm{gran}}=\mathrm{log}N_{2\mathrm{D}}\mathrm{log}N_{1\mathrm{D}}.$$
(11)
For simplicity, we used the standard HM model atmosphere (Holweger & Müller HM74 (1974); $`\xi _{\mathrm{mic}}=1`$ km/s) for this final step of the procedure ($`W_{1\mathrm{D}}N_{2\mathrm{D}}`$, $`W_{2\mathrm{D}}N_{1\mathrm{D}}`$). This is justified because all we need is a realistic curve-of-growth, $`W(N)`$, which is very similar for the averaged simulation snapshots and the HM atmosphere.

Case (A) and (B) give essentially the same corrections for weak lines, but the results may differ substantially for stronger lines, where $`W_{2\mathrm{D},\mathrm{A}}<W_{2\mathrm{D},\mathrm{B}}`$ (implying that the hydrodynamical velocity field provides an effective ‘microturbulence’ $`\xi _{\mathrm{mic}}(\mathrm{RHD})<1`$ km/s). For lines of arbitrary strength, we strongly favor the granulation corrections derived from case (B), because this strictly differential 1D/2D comparison is based on the same non-thermal velocity field, and thus cleanly separates the effect of the *temperature inhomogeneities* from that of the small-scale velocity field. In the following, $`\mathrm{\Delta }_{\mathrm{gran}}`$ therefore denotes the granulation abundance corrections derived from case (B).

For practical reasons, we have used only a small number of snapshots ($`S=11`$) for our present investigation. The statistical uncertainty of the derived abundance corrections is therefore appreciable. In order to quantify this uncertainty, we have also tabulated the standard deviation of the corrections $`\mathrm{\Delta }_{\mathrm{gran}}`$,

$$\sigma =\sqrt{\frac{1}{S1}\underset{s=1}{\overset{S}{}}\left(\mathrm{\Delta }_{\mathrm{gran},\mathrm{s}}\overline{\mathrm{\Delta }_{\mathrm{gran}}}\right)^2},$$
(12)
derived from the 11 individual snapshots. The quantity $`\sigma `$ can serve to judge the significance of the abundance corrections for each of the lines listed in the tables below.

## 5 Results for the Sun

We have applied the above procedure to a representative sample of lines, covering different elements, ionization stages, ionization energies, excitation energies, line strengths and wavelengths. In this differential study, we have adopted the same representative van der Waals broadening constant for all lines, $`\mathrm{log}C_6=30.3867`$ (definition according to Mihalas (Mi78 (1978)), which differs from that by Unsöld (U55 (1955)) by a factor $`2\pi `$), and have neglected line broadening due to electron collisions, $`C_4=0`$. The results are listed in Tables 1, 2, 3, and 4. Most of the lines are *’fictitious’ lines* (‘even’ wavelength entries in the tables), allowing a systematic study of the granulation effect as a function of the line parameters. The tables also contain a number of *’real’ lines* (table entries with ‘uneven’ wavelength) for Li, O, and Si, which have already been used with detailed abundance studies by Wedemeyer (SW01 (2001)) and Holweger (HH01 (2001)). In the latter work, the granulation abundance corrections were obtained by interpolation in the list of fictitious lines (Tables 1,3), which covers a sufficient range of the relevant line parameters.

### 5.1 Atoms with high ionization potential and ions

Atoms with high ionization potential like C i, N i, and O i, as well as singly ionized elements are representative of species in the main ionization stage. Neutral nitrogen (ionization potential $`\chi =14.5`$ eV) can be taken as the prototype of these species. Unlike C and O, N is not significantly affected by molecule formation or changes in the degree of ionization.

The corrections $`\mathrm{\Delta }_{\mathrm{gran}}`$ as a function of excitation energy, $`E_i`$, are displayed in Fig. 5 for our tabulated sample of weak fictitious lines ($`\lambda \mathrm{\hspace{0.17em}5500}`$ Å, $`W5`$ mÅ) of these species. Obviously, the corrections for high-excitation lines of the ‘majority species’ are similar for all elements, and can amount to -0.1 dex for weak lines with $`E_i10`$ eV. This conclusion is independent of the assumed velocity field which is irrelevant here: $`\mathrm{\Delta }_{\mathrm{gran},\mathrm{A}}\mathrm{\Delta }_{\mathrm{gran},\mathrm{B}}`$ for lines as weak as 5 mÅ (see tables).

#### 5.1.1 Dependence on line strength

The dependence of the granulation abundance corrections on line strength was investigated by test calculations for two additional sets of stronger N i lines. The two set of lines were obtained from the set of N i lines ($`\lambda \mathrm{\hspace{0.17em}5\hspace{0.17em}500}`$ Å) listed in Table 1 by increasing their equivalent widths, $`W_{1\mathrm{D}}`$, by a factor of 5 and 15, respectively. Thus, the test lines in first set of have $`W_{1\mathrm{D}}25`$ mÅ, those in the second set $`W_{1\mathrm{D}}80`$ mÅ. According to the results shown in Fig. 6, the granulation abundance corrections $`\mathrm{\Delta }_{\mathrm{gran}}`$ are systematically more positive for the stronger, otherwise identical lines. For $`E_i=12`$ eV the difference in $`\mathrm{\Delta }_{\mathrm{gran}}`$ between the lines with $`W_{1\mathrm{D}}80`$ mÅ and $`5`$ mÅ amounts to +0.07 dex.

#### 5.1.2 Dependence on wavelength

In order to get an idea about the dependence of the granulation abundance corrections on the wavelength of the spectral lines, we also performed test calculations for five additional sets of N i lines, obtained by shifting the weak N i lines listed in Table 1 from $`\lambda =\mathrm{5\hspace{0.17em}500}`$ Å to $`\lambda =\mathrm{4\hspace{0.17em}000}`$, $`\mathrm{8\hspace{0.17em}000}`$, $`\mathrm{10\hspace{0.17em}000}`$, $`\mathrm{12\hspace{0.17em}000}`$, and $`\mathrm{16\hspace{0.17em}000}`$ Å, respectively, enforcing the same *reduced* equivalent width, $`W/\lambda `$, i.e. $`W_{1\mathrm{D}}\lambda /\mathrm{5\hspace{0.17em}500}5`$ mÅ.

The results shown in Fig. 7 indicate that the granulation abundance corrections $`\mathrm{\Delta }_{\mathrm{gran}}`$ indeed depend on wavelength, in the sense that lines in the blue part of the spectrum show systematically more negative corrections than those in the red part. This general trend is evident over the whole wavelength range, apart from a slight reversal between $`\lambda \mathrm{\hspace{0.17em}8\hspace{0.17em}000}`$ and $`\mathrm{10\hspace{0.17em}000}`$ Å which is related to the Paschen jump at $`\lambda \mathrm{\hspace{0.17em}8208}`$ Å. For $`E_i=12`$ eV the difference in $`\mathrm{\Delta }_{\mathrm{gran}}`$ between a weak line at $`\lambda \mathrm{\hspace{0.17em}4\hspace{0.17em}000}`$ and $`\mathrm{10\hspace{0.17em}000}`$ Å amounts to -0.1 dex.

### 5.2 Atoms with low ionization potential

Elements with low ionization potential are predominantly singly ionized in the solar atmosphere. The particle density of the remaining neutral atoms of such elements is particularly sensitive to temperature variations. Molecules fall into the same category. Neutral iron (ionization potential $`\chi =7.9`$ eV) can be taken as representative for these species. According to our models, the fraction of neutral iron is nowhere exceeding 26%; it shows a maximum near $`\mathrm{log}\tau _{5500}=3`$ and has dropped to less than 3.5% at $`\mathrm{log}\tau _{5500}=0`$.

The corrections $`\mathrm{\Delta }_{\mathrm{gran}}`$ for a sample of weak fictitious lines ($`\lambda \mathrm{\hspace{0.17em}5500}`$ Å, $`W5`$ mÅ) of these ‘minority species’ are displayed in Fig. 8 as a function of excitation energy, $`E_i`$. We note that the largest granulation effects are found for low-excitation lines of neutral atoms with an ionization potential between 6 and 8 eV. The magnitude of the correction depends on the element and is as large as -0.3 dex for Ti i, $`E_i=0`$ eV. Again, we point out that this result does not depend on the assumed non-thermal velocity field.

It is interesting to note that Asplund (2000a ), using 3D models, finds corrections in the opposite direction for the ’real’ Si I lines listed in Table 2. We emphasize that this is not an indication of fundamentally different “granulation corrections” in 2D and 3D. The corrections given by Asplund with respect to the Holweger-Müller (HM74 (1974)) model include not only the effect of temperature inhomogeneities but also the influence of the different mean stratifications. Obviously, the difference in the mean structure is the more important of the two opposing effects in this case. Indeed, we find from our own models a *total* correction of -0.05 dex with respect to the Holweger-Müller model, compared to a the tabulated ”granulation correction” of +0.02 dex, explaining the apparently inconsistent results.

#### 5.2.1 Dependence on line strength

As for N i, the dependence of the granulation abundance corrections on line strength was investigated by test calculations for two additional sets of stronger Fe i lines, obtained from the set of Fe i lines listed in Table 3 by increasing their equivalent widths, $`W_{1\mathrm{D}}`$, by a factor of 5 and 10, respectively. The test lines in first set of have then $`W_{1\mathrm{D}}35`$ mÅ, those in the second set $`W_{1\mathrm{D}}70`$ mÅ. The results are shown in Fig. 9. Qualitatively, the granulation abundance corrections $`\mathrm{\Delta }_{\mathrm{gran}}`$ for Fe i show the same dependence on line strength as already seen for N i: they are systematically more positive for the stronger, otherwise identical lines.

The corrections for $`E_i=0`$ and $`1`$ eV, $`W_{1\mathrm{D}}70`$ mÅ, shown in Fig. 9 are certainly too small (too positive), because a significant part of the line absorption comes from layers located outside the model atmosphere on which the line formation calculations are based ($`\tau _{5500}<10^5`$). This situation results in a stronger underestimation of the equivalent width for the inhomogeneous case. Hence, $`\mathrm{\Delta }_{\mathrm{gran}}`$ should actually be more negative, and we can expect abundance corrections $`\mathrm{\Delta }_{\mathrm{gran}}<0.1`$ dex even for the stronger ground state lines.

#### 5.2.2 Dependence on wavelength

We have performed test calculations for five additional sets of Fe i lines, obtained by shifting the standard Fe i lines listed in Table 3 from $`\lambda =\mathrm{5\hspace{0.17em}500}`$ Å to $`\lambda =\mathrm{4\hspace{0.17em}000}`$, $`\mathrm{8\hspace{0.17em}000}`$, $`\mathrm{10\hspace{0.17em}000}`$, $`\mathrm{12\hspace{0.17em}000}`$, and $`\mathrm{16\hspace{0.17em}000}`$ Å, respectively, enforcing the same *reduced* equivalent width, $`W/\lambda `$, i.e. $`W_{1\mathrm{D}}\lambda /\mathrm{5\hspace{0.17em}500}5`$ mÅ.

The results shown in Fig. 10 indicate that the wavelength dependence of the granulation abundance corrections $`\mathrm{\Delta }_{\mathrm{gran}}`$ for Fe i has the opposite sign as in the case of N i. Lines in the blue part of the spectrum tend to show smaller (less negative) corrections than those at longer wavelengths. This general trend is clearly evident from $`\lambda \mathrm{\hspace{0.17em}4000}`$ to $`8000`$ Å; at longer wavelengths the variation of $`\mathrm{\Delta }_{\mathrm{gran}}`$ with $`\lambda `$ is much reduced and no longer strictly monotonic (for an explanation see Paper II).

We have to point out that the results for $`\lambda `$ 16 000 Å are uncertain, because a significant part of the line absorption comes from layers above $`\tau _{16000}<10^5`$, so the line formation region is not entirely covered by our model atmosphere.

## 6 Conclusions

In an effort to estimate the - hitherto largely unknown - effects of photospheric temperature fluctuations on spectroscopic abundance determinations, we have carried out numerical simulations of the problem of *LTE line formation in the inhomogeneous solar photosphere* based on detailed 2-dimensional radiation hydrodynamics models of the convective surface layers of the Sun.

For a variety of spectral lines of different elements, we have computed synthetic line profiles from the inhomogeneous 2D hydrodynamical atmosphere and from the corresponding 1D plane-parallel model, respectively. By means of a strictly differential 1D/2D comparison of the emergent equivalent widths, we have derived so called “granulation abundance corrections” for the individual lines, which have to be applied to standard abundance determinations based on homogeneous 1D model atmospheres in order to correct for the influence of the photospheric temperature fluctuations. The ‘classical’ problem of finding the most appropriate mean vertical temperature stratification is *not* addressed here.

Using the concept of “fictitious” spectral lines, we were able to investigate systematically the dependence of the “granulation abundance corrections” on the basic line parameters like excitation potential, line strength, wavelength, element and ionization stage. For many lines of practical relevance, it should be possible to estimate the magnitude of the abundance correction by interpolation in the graphs and tables provided in this paper. This approach may often be an acceptable alternative to a detailed fitting of individual line profiles based on hydrodynamical simulations.

In general, we find a *line strengthening* in the presence of temperature inhomogeneities, implying mostly *negative* “granulation abundance corrections”, i.e. standard analysis based on plane-parallel atmospheres tends to *overestimate* abundances. The physical reason for the line strengthening is primarily the non-linear temperature dependence of the line opacity due to thermal ionization and excitation, as demonstrated in Paper II.

One remarkable result is that all lines of our sample with an excitation potential around $`E_i5`$ eV are practically *insensitive* to granulation effects, regardless of element and ionization stage. Moderate granulation corrections ($`\mathrm{\Delta }_{\mathrm{gran}}0.1`$ dex) are found for weak, high-excitation lines ($`E_i10`$ eV) of ions and atoms with high ionization potential like N i. The largest corrections are found for ground state lines ($`E_i=0`$ eV) of neutral atoms with an ionization potential between 6 and 8 eV like Mg i, Ca i, Ti i, Fe i, amounting to $`\mathrm{\Delta }_{\mathrm{gran}}0.3`$ dex in the case of Ti i.

For given excitation potential, the granulation corrections are systematically more positive for stronger lines. The wavelength dependence of $`\mathrm{\Delta }_{\mathrm{gran}}`$, however, depends on the type of line (‘majority’ or ‘minority’ species).

We recall that the corrections $`\mathrm{\Delta }_{\mathrm{gran}}`$ account only for the influence of the temperature and pressure fluctuations relative to a *given mean structure*. We have suppressed possible effects related to 1D/2D differences in the details of the non-thermal velocity field. Such additional effects, which apply to partly saturated lines only, should be small, however, since they are largely eliminated by the proper choice of the microturbulence parameter $`\xi _{\mathrm{mic}}`$.

We point out that the magnitude of the “granulation abundance corrections” derived in this work should be considered as upper limits ($`|\mathrm{\Delta }_{\mathrm{gran}}|(2D,LTE)>|\mathrm{\Delta }_{\mathrm{gran}}|(3D,NLTE)`$), because (i) possible NLTE-effects, which were ignored in this study, tend to reduce the fluctuations of the line opacity (e.g. Kiselman Ki97 (1997), Cayrel & Steffen CS00 (2000), Asplund 2000b ), and (ii) the amplitude of the temperature fluctuations is systematically overestimated in 2D relative to 3D convection models (see Fig. 4). The large negative corrections found for the low-excitation lines of atoms with low ionization potential are the result of the strongly increasing amplitude of the temperature fluctuations towards the upper photosphere as predicted by our present 2D hydrodynamical simulations. In a more realistic 3D atmosphere, this feature will therefore be less pronounced.

Even though the results are still somewhat uncertain, the present study should be helpful in providing upper bounds for possible errors of spectroscopic abundance analyses, and for selecting spectral lines which are as insensitive as possible to the effects of photospheric temperature inhomogeneities.

For a physical explanation of the numerical results obtained in this work, a detailed investigation of the process of LTE line formation in an inhomogeneous stellar atmosphere based on the analysis of the transfer equation and the evaluation of line depression contribution functions is presented in the second paper of this series (Paper II). With the help of this formalism, we can actually tell how the different atmospheric layers contribute to the granulation correction for any particular spectral line.

In the near future, we plan to extend the computation of abundance corrections to full 3D simulations for the Sun and other types of stars where convection extends into higher photospheric layers and the spectroscopic granulation effects may well be even larger.

###### Acknowledgements.

The 2D numerical convection simulations and line formation calculations were carried out on the CRAY T94 and CRAY SV1, respectively, at the *Rechenzentrum der Universität Kiel*. We thank the referee, Paul Barklem, for constructive criticism which helped to significantly improve the contents and presentation of this work.

## 7 Granulation Abundance Correction Tables