text:
# Energy-momentum Tensors in Matrix Theory and in Noncommutative Gauge Theories

## 1. Introduction

Noncommutative gauge theories can be realized by considering branes in string theory with a strong NS-NS two-form field . In our previous paper , we derived the energy-momentum tensors of these theories by computing disk amplitudes with one closed string and an arbitrary number of open strings and by taking the Seiberg-Witten limit  of these amplitudes. We found that the energy-momentum tensors involve the open Wilson lines .111Various aspects of the open Wilson lines and related issues were discussed in . However, they do not reduce to the ones in the commutative theories in the limit $`\theta 0`$, where $`\theta `$ is the noncommutative parameter. This is because the Seiberg-Witten limit does not commute with the commutative limit.222In this paper, we study the energy-momentum tensor coupled to gravity in the bulk. The energy-momentum tensor derived in the Nöther procedure, such as the one discussed in , would couple to the open string metric on the branes. We also found that the energy-momentum tensors are conserved in interesting ways. In particular, in theories derived from bosonic string, the energy-momentum tensors are kinematically conserved, namely, the conservation law holds identically for any field configuration irrespective of the equations of motion.

It turns out that we can use the same method to derive the energy-momentum tensor of Matrix Theory .333See also a recent review . If we perform the dimensional reduction of the noncommutative theories along their noncommutative directions, the Seiberg-Witten limit reduces to the DKPS limit  (or the Sen-Seiberg limit  in the context of Matrix Theory ). Therefore, the energy-momentum tensor of Matrix Theory can be derived as a special case of the computation in  where there are no noncommutative directions. In , the energy-momentum tensor of Matrix Theory was deduced from computations of one-loop amplitudes in Matrix Theory and from their comparison with graviton exchange amplitudes. The energy-momentum tensor we derive from the disk amplitudes perfectly agrees with that obtained in  including the structure of the higher moments.

It has been shown in  and  that noncommutative gauge theories are realized as certain backgrounds of Matrix Theory and of the IIB matrix model. This implies a natural correspondence of gauge invariant observables between Matrix Theory and noncommutative gauge theories. In fact, it was in this context  that the open Wilson lines were discovered as a basis for observables in noncommutative theories. In , it was proposed to use this correspondence to derive the energy-momentum tensors of the noncommutative theories.444See also . In this paper, we carry out this proposal explicitly. We find that the resulting energy-momentum tensors agree exactly with those derived in , both in the superstring case and in the bosonic string case. The correspondence between Matrix Theory and noncommutative theories also leads to a simpler proof of the conservation law of the energy-momentum tensors, which came out rather miraculously in our earlier work .

This paper is organized as follows. In Section 2, we will derive the energy-momentum tensor of Matrix Theory by computing disk amplitudes with one closed string and an arbitrary number of open strings and by taking the DKPS limit. As we mentioned in the above, this can be regarded as a special case of the computation in . In Section 3, we use the correspondence between Matrix Theory and noncommutative gauge theories to derive the energy-momentum tensors in noncommutative theories and find an agreement with our earlier result. We will study the case of bosonic string in Section 4.

Let us summarize the notations used in this paper. The coordinates in the bulk are denoted by $`x^M`$. For the discussion regarding the DKPS limit of D-branes, we use the late Greeks ($`\mu ,\nu ,\mathrm{}`$) for the directions along the branes and $`I,J,\mathrm{}`$ for those transverse to the branes:

(1)
$$x^M=\{\begin{array}{cc}x^\mu \hfill & \text{for the directions along the branes},\hfill \\ x^I\hfill & \text{for the transverse directions}.\hfill \end{array}$$

In the context of noncommutative gauge theories, we use the romans ($`i,j,\mathrm{}`$) for the noncommutative directions on the branes, the late Greeks ($`\mu ,\nu ,\mathrm{}`$) for the commutative directions, and the early Greeks ($`\alpha ,\beta ,\mathrm{}`$) for the directions transverse to the branes:

(2)
$$x^M=\{\begin{array}{cc}x^\mu \hfill & \text{for the commutative directions along the branes},\hfill \\ x^i\hfill & \text{for the noncommutative directions along the branes},\hfill \\ x^\alpha \hfill & \text{for the transverse directions}.\hfill \end{array}$$

As we will see, the transverse directions $`x^I`$ transmute into $`x^i`$ and $`x^\alpha `$.

In the noncommutative directions, the coordinates obey the Heisenberg relation,

(3)
$$[x^i,x^j]=i\theta ^{ij}.$$

The closed string metric is denoted by $`g_{MN}`$ and the open string metric $`G_{MN}`$ in the zero-slope limit is given by

(4)
$$G^{ij}=\frac{1}{(2\pi \alpha ^{})^2}\theta ^{im}\theta ^{jn}g_{mn},G^{\mu \nu }=g^{\mu \nu }.$$

The noncommutative gauge theory is defined by taking the zero-slope limit, $`\alpha ^{}0`$, while keeping $`G^{ij}`$ and $`\theta ^{ij}`$ finite (the Seiberg-Witten limit ). This means that the closed string metric along the noncommutative direction $`g_{ij}`$ is scaled as $`\alpha _{}^{}{}_{}{}^{2}`$.

## 2. Energy-momentum tensor of Matrix Theory

In , the energy-momentum tensors of noncommutative gauge theories are derived by computing the disk amplitudes with one closed string and an arbitrary number of open strings. It turns out that this method is applicable to the computation of the actions of multiple D-branes in weakly curved background, or equivalently, the energy-momentum tensor of Matrix Theory.

The action of a single D$`p`$-brane in curved space is given by

(5)
$$S=\frac{1}{g_{YM}^2}𝑑x\sqrt{detg}\left[\frac{1}{4}g^{\mu \nu }g^{\rho \sigma }F_{\mu \rho }F_{\nu \sigma }+\frac{1}{2(2\pi \alpha ^{})^2}g^{\mu \nu }g_{IJ}_\mu Y^I_\nu Y^J\right],$$

where the scalar fields $`Y^I(x)`$ describe the location of the brane in the transverse directions and the Yang-Mills coupling constant $`g_{YM}`$ is given in terms of the string coupling constant $`g_s`$ and $`\alpha ^{}`$ as follows:

(6)
$$g_{YM}^2=g_s(2\pi )^{p2}\alpha _{}^{}{}_{}{}^{\frac{p3}{2}}.$$

The integral $`𝑑x`$ is over the $`p+1`$ dimensional worldvolume of the D-brane. The action (5) is obtained as a zero-slope limit $`\alpha ^{}0`$ known as the DKPS limit  (or the Sen-Seiberg limit  in the context of Matrix Theory ) given by

(7)
$$Y^IO(\alpha ^{}),k_IO(\alpha _{}^{}{}_{}{}^{1}),$$

while the other quantities $`A_\mu ,k_\mu ,g_{\mu \nu },g_{IJ}`$ and $`g_{YM}`$ are kept finite as $`\alpha ^{}0`$. Here $`k_M=(k_\mu ,k_I)`$ is the Fourier mode conjugate to the coordinate $`x^M=(x^\mu ,x^I)`$. In this limit, $`\mathrm{\Phi }^I`$ defined by

(8)
$$\mathrm{\Phi }^I(x)=\frac{Y^I(x)}{2\pi \alpha ^{}},$$

is kept finite. In the case of multiple D-branes, where $`A_\mu `$ and $`Y^I`$ become matrix valued, a complete action in curved space is not known although its various aspects have been studied in  and . Here we will study the first term in the expansion,

(9)
$$g_{IJ}(x,Y(x))=\eta _{IJ}+h_{IJ}\mathrm{exp}\left[ik_\mu x^\mu +ik_IY^I(x)\right]+O(h^2).$$

The DKPS limit (7) defined in this way is the same as the Seiberg-Witten limit of noncommutative gauge theories as far as the commutative directions are concerned. Therefore, we can derive the energy-momentum tensors of multiple D-branes by considering the energy-momentum tensors of noncommutative gauge theories found in  and by taking their dimensional reductions in the noncommutative directions.

It would be instructive to recapitulate the derivation of the energy-momentum tensor in , in the language of Matrix Theory. The worldsheet action is

(10)
$$S=\frac{1}{2\pi }d^2\sigma \left[\frac{2}{\alpha ^{}}g_{MN}X^M\overline{}X^N+\frac{1}{\alpha ^{}}g_{MN}\psi ^M\overline{}\psi ^N+\frac{1}{\alpha ^{}}g_{MN}\overline{\psi }^M\overline{\psi }^N\right].$$

When the worldsheet disk is represented by the upper-half plane of the complex plane, the propagators for $`X^M(z)`$ are

(11) $`X^M(z)X^N(w)`$
$`=\alpha ^{}\left[g^{MN}\mathrm{log}|zw|g^{MN}\mathrm{log}|z\overline{w}|+2G^{MN}\mathrm{log}|z\overline{w}|\right],`$
and those for the fermions are

(12) $`\psi ^M(z)\psi ^N(w)`$ $`={\displaystyle \frac{\alpha ^{}}{zw}}g^{MN},`$
$`\psi ^M(z)\overline{\psi }^N(\overline{w})`$ $`={\displaystyle \frac{\alpha ^{}}{z\overline{w}}}\left(g^{MN}+2G^{MN}\right),`$
$`\overline{\psi }^M(\overline{z})\psi ^N(w)`$ $`={\displaystyle \frac{\alpha ^{}}{\overline{z}w}}\left(g^{MN}+2G^{MN}\right),`$
$`\overline{\psi }^M(\overline{z})\overline{\psi }^N(\overline{w})`$ $`={\displaystyle \frac{\alpha ^{}}{\overline{z}\overline{w}}}g^{MN},`$
where

(13)
$$G^{\mu \nu }=g^{\mu \nu },G^{IJ}=0.$$

The fermions on the boundary $`\mathrm{\Psi }^M(t)`$ are defined by

(14)
$$\mathrm{\Psi }^\mu (t)=\frac{1}{2}\left[\psi ^\mu (t)+\overline{\psi }^\mu (t)\right],\mathrm{\Psi }^I=\frac{1}{2}\left[\psi ^I(t)\overline{\psi }^I(t)\right],$$

and the boundary-boundary and bulk-boundary propagators for the fermions are given by

(15) $`\mathrm{\Psi }^\mu (t)\mathrm{\Psi }^\nu (t^{})`$ $`={\displaystyle \frac{\alpha ^{}}{tt^{}}}g^{\mu \nu },`$ $`\mathrm{\Psi }^I(t)\mathrm{\Psi }^J(t^{})`$ $`={\displaystyle \frac{\alpha ^{}}{tt^{}}}g^{IJ},`$
(16) $`\psi ^\mu (z)\mathrm{\Psi }^\nu (t)`$ $`={\displaystyle \frac{\alpha ^{}}{zt}}g^{\mu \nu },`$ $`\overline{\psi }^\mu (\overline{z})\mathrm{\Psi }^\nu (t)`$ $`={\displaystyle \frac{\alpha ^{}}{\overline{z}t}}g^{\mu \nu },`$
$`\psi ^I(z)\mathrm{\Psi }^J(t)`$ $`={\displaystyle \frac{\alpha ^{}}{zt}}g^{IJ},`$ $`\overline{\psi }^I(\overline{z})\mathrm{\Psi }^J(t)`$ $`={\displaystyle \frac{\alpha ^{}}{\overline{z}t}}g^{IJ},`$
where $`t`$ and $`t^{}`$ are on the boundary.

It is convenient to use the vertex operator in the $`(1,1)`$-picture for the closed string,

(17)
$$V^{(1,1)}(z)=\frac{1}{2}\delta (\gamma )\delta (\overline{\gamma })h_{MN}(k)\psi ^M(z)\overline{\psi }^N(\overline{z})e^{ikX(z)},$$

where $`\gamma `$ and $`\overline{\gamma }`$ are bosonic ghosts, and the operator in the $`0`$-picture for open string,

(18)
$$\begin{array}{c}U^{(0)}(t)=A_\mu (X)\frac{dX^\mu }{dt}+\mathrm{\Phi }^I(X)ig_{IJ}_{}X^JF_{\mu \nu }(X)\mathrm{\Psi }^\mu \mathrm{\Psi }^\nu \hfill \\ \hfill 2D_\mu \mathrm{\Phi }^I(X)g_{IJ}\mathrm{\Psi }^\mu \mathrm{\Psi }^Jig_{IK}g_{JL}[\mathrm{\Phi }^I(X),\mathrm{\Phi }^J(X)]\mathrm{\Psi }^K\mathrm{\Psi }^L.\end{array}$$

Here we used $`\mathrm{\Phi }^I`$ defined in (8) instead of $`Y^I`$ to make it clear that the open string vertex operator is $`O(1)`$ in the zero-slope limit. The correlation function which we need to evaluate is then

(19)
$$(z\overline{z})^2V^{(1,1)}(z)trP\mathrm{exp}\left(i_{\mathrm{}}^{\mathrm{}}𝑑tU^{(0)}(t)\right),$$

with a proper gauge-fixing of the $`SL(2,R)`$ invariance. Here the symbol $`P`$ denotes the path-ordering with respect to $`t`$.

By definition, all the propagators are $`O(\alpha ^{})`$ whereas the open vertex operator is $`O(1)`$ in the zero-slope limit. Therefore, it costs $`\alpha ^{}`$ each time we use a propagator. However, when we make a contraction between an open string vertex operator and the $`e^{ik_IX^I}`$ part of the graviton vertex operator, the factor $`\alpha ^{}`$ is cancelled by the fact that the momentum $`k_I`$ scales as $`O(\alpha _{}^{}{}_{}{}^{1})`$, so the contraction remains finite for any number of open string vertex operators.

(20) $`e^{ikX(z)}tr{\displaystyle \underset{a}{}}\mathrm{\Phi }^I(X(t_a))ig_{IJ}_{}X^J(t_a)`$
$`={\displaystyle 𝑑xe^{ik_\mu x^\mu }tr\underset{a}{}2\pi \alpha ^{}k_I\mathrm{\Phi }^I(x)\frac{\tau (t_a,z)}{t_a}}.`$
Here the function $`\tau (t,z)`$ is given by

(21)
$$\tau (t,z)=\frac{1}{2\pi i}\mathrm{log}\left(\frac{tz}{t\overline{z}}\right),$$

and we used that

(22)
$$i_{}\mathrm{log}|zt|^2=\frac{z\overline{z}}{(zt)(\overline{z}t)}=2\pi i\frac{\tau (t,z)}{t}.$$

For a fixed value of $`z`$, $`\tau (t,z)`$ is a monotonically increasing function of $`t`$ and

(23)
$$\tau (\mathrm{},z)\tau (\mathrm{},z)=1.$$

Therefore, we have

(24) $`e^{ikX(z)}trP\mathrm{exp}\left(i{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t\mathrm{\Phi }^I(X(t))ig_{IJ}_{}X^J(t)\right)`$
$`={\displaystyle 𝑑xe^{ik_\mu x^\mu }trP\mathrm{exp}\left(i_0^1𝑑\tau 2\pi \alpha ^{}k_I\mathrm{\Phi }^I(x)\right)}.`$
Note that the path-ordering with respect to $`\tau `$ on the right-hand side is inherited from that with respect to $`t`$ on the left-hand side because $`\tau `$ is a monotonically increasing function of $`t`$. Note also that the $`A_\mu (X)dX^\mu /dt`$ part of the open string vertex operator (18) does not contribute to the energy-momentum tensor in the zero-slope limit because there is no divergent factor like $`k_I`$ which can compensate the factor of $`\alpha ^{}`$ coming from the propagator.

We also need to contract the fermions in the closed string vertex operator to obtain a non-vanishing result. The contribution to the graviton vertex operator starts at $`O(\alpha _{}^{}{}_{}{}^{3})`$ with three contractions of the fermions.

(25) $`(z\overline{z}){\displaystyle \frac{1}{2}}\left[\psi ^I(z)\overline{\psi }^J(\overline{z})+\psi ^J(z)\overline{\psi }^I(\overline{z})\right]trP\mathrm{exp}\left(i{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑tU^{(0)}(t)\right)`$
$`=8\pi ^2\alpha _{}^{}{}_{}{}^{3}tr[g^{\mu \nu }{\displaystyle _0^1}d\tau _1D_\mu \mathrm{\Phi }^I(x){\displaystyle _0^1}d\tau _2D_\nu \mathrm{\Phi }^J(x)`$
$`+g_{KL}{\displaystyle _0^1}d\tau _1[\mathrm{\Phi }^I(x),\mathrm{\Phi }^K(x)]{\displaystyle _0^1}d\tau _2[\mathrm{\Phi }^L(x),\mathrm{\Phi }^J(x)]]+O(\alpha _{}^{}{}_{}{}^{4})`$
$`+\text{the dilaton part}.`$
After taking into account the contributions from the ghosts, we obtain the components $`T^{IJ}`$ of the energy-momentum tensor in the case of superstring as follows:

(26)
$$\begin{array}{c}2\alpha ^{}dxe^{ikx}tr[\mathrm{exp}\left(i_0^1d\tau k_IY^I(x)\right)\hfill \\ \hfill \times [g^{\mu \nu }_0^1d\tau _1D_\mu Y^I(x)_0^1d\tau _2D_\nu Y^J(x)\\ \hfill +\frac{g_{KL}}{(2\pi \alpha ^{})^2}_0^1d\tau _1[Y^I(x),Y^K(x)]_0^1d\tau _2[Y^L(x),Y^J(x)]]],\end{array}$$

where we used $`Y^I`$ instead of $`\mathrm{\Phi }^I`$ in (8). Note that the ordering of the fields is determined by the values of $`\tau `$, $`\tau _1`$ and $`\tau _2`$ which are inherited from the positions of open string vertex operators. We can transform the expression for $`T^{IJ}`$ into a more convenient form. For any operators $`A`$ and $`B`$,

(27) $`tr\left[\mathrm{exp}\left(i{\displaystyle _0^1}𝑑\tau kY\right){\displaystyle _0^1}𝑑\tau _1A{\displaystyle _0^1}𝑑\tau _2B\right]`$
$`={\displaystyle _0^1}𝑑\tau _1{\displaystyle _{\tau _1}^1}𝑑\tau _2trAe^{i(\tau _2\tau _1)kY}Be^{i(1+\tau _1\tau _2)kY}`$
$`+{\displaystyle _0^1}𝑑\tau _1{\displaystyle _0^{\tau _1}}𝑑\tau _2trAe^{i(1\tau _1+\tau _2)kY}Be^{i(\tau _1\tau _2)kY}`$
$`={\displaystyle _0^1}𝑑\tau \left[{\displaystyle _0^{1\tau }}𝑑\tau ^{}+{\displaystyle _{1\tau }^1}𝑑\tau ^{}\right]trAe^{i\tau kY}Be^{i(1\tau )kY}`$
$`={\displaystyle _0^1}𝑑\tau trAe^{i\tau kY}Be^{i(1\tau )kY}.`$
Therefore, we can write

(28)
$$\begin{array}{c}T^{IJ}=2\alpha ^{}dxe^{ikx}tr_0^1d\tau [g^{\mu \nu }D_\mu Y^Ie^{i\tau k_IY^I}D_\nu Y^Je^{i(1\tau )k_IY^I}\hfill \\ \hfill +\frac{g_{KL}}{(2\pi \alpha ^{})^2}[Y^I,Y^K]e^{i\tau k_IY^I}[Y^L,Y^J]e^{i(1\tau )k_IY^I}].\end{array}$$

Furthermore, we can carry out the remaining $`\tau `$-integral and the result is expressed in terms of the symmetrized trace .555 We thank S. Das, W. Taylor and S. Trivedi for discussion on this point. To see this we note that, for any operators $`A`$ and $`B`$,

(29) $`{\displaystyle _0^1}𝑑\tau tr\left[Ae^{i\tau kY}Be^{i(1\tau )kY}\right]`$
$`={\displaystyle \underset{n,m=0}{\overset{\mathrm{}}{}}}{\displaystyle _0^1}𝑑\tau tr\left[A{\displaystyle \frac{1}{n!}}(i\tau kY)^nB{\displaystyle \frac{1}{m!}}(i(1\tau )kY)^m\right]`$
$`={\displaystyle \underset{n,m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{n!m!}{(n+m+1)!}}tr\left[A{\displaystyle \frac{1}{n!}}(ikY)^nB{\displaystyle \frac{1}{m!}}(ikY)^m\right]`$
$`={\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{p!}}{\displaystyle \frac{1}{p+1}}{\displaystyle \underset{q=0}{\overset{p}{}}}tr\left[A(ikY)^qB(ikY)^{pq}\right].`$
Since

(30) $`Str\left(C^pAB\right)`$ $`={\displaystyle \frac{1}{(p+1)!}}{\displaystyle \underset{q=0}{\overset{p}{}}}p!tr\left[AC^qBC^{pq}\right]`$
$`={\displaystyle \frac{1}{p+1}}{\displaystyle \underset{q=0}{\overset{p}{}}}tr\left[AC^qBC^{pq}\right],`$
where $`Str`$ is the symmetrized trace, we can write

(31)
$$_0^1𝑑\tau tr\left[Ae^{i\tau kY}Be^{i(1\tau )kY}\right]=Str\left(e^{ikY}AB\right).$$

Thus $`T^{IJ}`$ is expressed as

(32)
$$\begin{array}{c}T^{IJ}=2\alpha ^{}dxe^{ikx}Str[e^{ik_IY^I}g^{\mu \nu }D_\mu Y^ID_\nu Y^J\hfill \\ \hfill +e^{ik_IY^I}\frac{g_{KL}}{(2\pi \alpha ^{})^2}[Y^I,Y^K][Y^L,Y^J]],\end{array}$$

where we symmetrized over all possible orderings of $`Y^I`$, $`D_\mu Y^I`$ and $`[Y^I,Y^J]`$ inside the trace. Similarly, we can derive the expressions for $`T^{\mu J}`$ and $`T^{\mu \nu }`$.

(33) $`T^{\mu J}`$ $`=2\alpha ^{}g^{\mu \mu ^{}}{\displaystyle 𝑑xe^{ikx}tr\left[e^{ik_IY^I}D_\mu ^{}Y^J\right]},`$
(34) $`T^{\mu \nu }`$ $`=2\alpha ^{}g^{\mu \nu }{\displaystyle 𝑑xe^{ikx}tr\left[e^{ik_IY^I}\right]}.`$
These expressions agree with those derived in  using one-loop amplitudes of Matrix Theory. Our result for the graviton can be extended to all other closed string states including massive ones as discussed in .

The appearance of the symmetrized trace is not an obvious consequence of string theory computation. To the contrary, in the case of bosonic string, the ordering of operators in the energy-momentum tensor does not obey the symmetrized trace prescription as we will see in Section 4.

The conservation of the energy-momentum tensor,

(35)
$$k_\mu T^{\mu J}+k_IT^{IJ}=0,$$

has been verified by Van Raamsdonk . For the D$`(1)`$-branes, it goes as

(36) $`k_IT^{IJ}`$ $`tr{\displaystyle _0^1}𝑑\tau [kY,Y^K]e^{i\tau kY}[Y_K,Y^J]e^{i(1\tau )kY}`$
$`=itr{\displaystyle _0^1}𝑑\tau {\displaystyle \frac{d}{d\tau }}\left(Y^Ke^{i\tau kY}[Y_K,Y^J]e^{i(1\tau )kY}\right)`$
$`=itr[Y^K,[Y_K,Y^J]]e^{ikY}`$
$`=0,`$
where $`kYk_IY^I`$ and $`Y_Kg_{KL}Y^L`$. In the last line, we used the equation of motion of Matrix Theory,

(37)
$$[Y^K,[Y_K,Y^J]]=0.$$

It is easily generalized to the case of D$`p`$-branes using the identity

(38)
$$D_\mu (e^A)=_0^1𝑑\tau e^{(1\tau )A}D_\mu Ae^{\tau A}.$$

## 3. Relation to noncommutative gauge theory

Given the energy-momentum tensor of Matrix Theory, we can derive those of noncommutative gauge theories following the suggestion in . Let us carry this out explicitly here. We divide the transverse directions $`x^I`$ into two, $`x^I=(x^\alpha ,x^i)`$. To fit with the standard notation in noncommutative gauge theories, we rescale the coordinates $`x^i`$ so that $`g_{ij}O(\alpha _{}^{}{}_{}{}^{2})`$. On the other hand, $`g_{\alpha \beta }`$ as well as $`g_{\mu \nu }`$ remains finite. Consider the background

(39)
$$Y^\alpha =0,Y^i=x^i\text{such that}[x^i,x^j]=i\theta ^{ij},$$

and expand the scalar fields around this background as follows:

(40)
$$Y^\alpha =2\pi \alpha ^{}\mathrm{\Phi }^\alpha O(\alpha ^{}),Y^i=x^i+\theta ^{ij}A_jO(1).$$

We take the limit $`k_\alpha O(\alpha _{}^{}{}_{}{}^{1})`$ and $`g_{ij}O(\alpha _{}^{}{}_{}{}^{2})`$ keeping $`(\mathrm{\Phi }^\alpha ,A^i,A_\mu )`$, $`(k_i,k_\mu )`$ and $`(g_{\alpha \beta },G_{\mu \nu }=g_{\mu \nu })`$ finite. This is identical to the limit taken in  following . An open Wilson line, which constitutes a basic ingredient of the energy-momentum tensors of the noncommutative gauge theories, emerges in this background as follows :

(41)
$$\begin{array}{c}tr\mathrm{exp}\left[ik_IY^I\right]=tr\mathrm{exp}\left[ik_ix^i+il^iA_i+iy_\alpha \mathrm{\Phi }^\alpha \right]\hfill \\ \hfill 𝑑x\left[\mathrm{exp}\left(i_0^1𝑑\tau \left(l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau )\right)\right)e^{ik_ix^i}\right],\end{array}$$

where $`l^i=k_j\theta ^{ji}`$, $`y_\alpha =2\pi \alpha ^{}k_\alpha `$. Here we mapped the matrices into the functions on noncommutative space with the star product on the right-hand side. The symbol $`[\mathrm{}]`$ means that we take the star product in the expression in $`[\mathrm{}]`$ taking into account the path-ordering along the open Wilson line between $`x^i`$ and $`x^i+l^i`$. The integral $`𝑑x`$ here is over the noncommutative directions $`x^i`$ along the branes which constitutes the integral over the whole worldvolume of the branes together with the previous integral over the commutative directions $`x^\mu `$. In the zero-slope limit, the on-shell condition of the graviton is

(42)
$$G_{ij}l^il^j+g^{\alpha \beta }y_\alpha y_\beta =0.$$

For the cases with operator insertions along the Wilson line, the following formulas are useful:

(43) $`D_\mu Y^I`$ $`=\{\begin{array}{cc}2\pi \alpha ^{}D_\mu \mathrm{\Phi }^\alpha \hfill & \text{when }I=\alpha ,\hfill \\ \theta ^{ij}F_{\mu j}\hfill & \text{when }I=i,\hfill \end{array}`$
(44) $`[Y^I,Y^J]`$ $`=\{\begin{array}{cc}(2\pi \alpha ^{})^2[\mathrm{\Phi }^\alpha ,\mathrm{\Phi }^\beta ]\hfill & \text{when }(I,J)=(\alpha ,\beta ),\hfill \\ 2\pi i\alpha ^{}\theta ^{jj^{}}D_j^{}\mathrm{\Phi }^\alpha \hfill & \text{when }(I,J)=(\alpha ,j),\hfill \\ i\theta ^{ii^{}}\theta ^{jj^{}}(F_{i^{}j^{}}\theta _{i^{}j^{}}^1)\hfill & \text{when }(I,J)=(i,j),\hfill \end{array}`$
and

(45)
$$\begin{array}{c}tr_0^1𝑑\tau Ae^{i\tau k_IY^I}Be^{i(1\tau )k_IY^I}\hfill \\ \hfill dx[\mathrm{exp}\left(i_0^1d\tau (l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\right)\\ \hfill \times _0^1d\tau _1A(x+l\tau _1)_0^1d\tau _2B(x+l\tau _2)e^{ik_ix^i}].\end{array}$$

It is straightforward to verify that all the components of energy-momentum tensor of the noncommutative gauge theory in the case of superstring derived in  are correctly reproduced including the relative coefficients. For example, the expression for $`T^{ij}`$ is given by

(46)
$$\begin{array}{c}T^{ij}=2\alpha ^{}\theta ^{ii^{}}\theta ^{jj^{}}\hfill \\ \hfill \times dx[e^{ikx}\mathrm{exp}\left(i_0^1d\tau (l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\right)\\ \hfill \times _0^1d\tau _1_0^1d\tau _2[G^{mn}(F_{i^{}m}(x+l\tau _1)\theta _{i^{}m}^1)(F_{j^{}n}(x+l\tau _2)\theta _{j^{}n}^1)\\ \hfill +G^{\mu \nu }F_{i^{}\mu }(x+l\tau _1)F_{j^{}\nu }(x+l\tau _2)+g_{\alpha \beta }D_i^{}\mathrm{\Phi }^\alpha (x+l\tau _1)D_j^{}\mathrm{\Phi }^\beta (x+l\tau _2)]],\end{array}$$

and the other components are presented in the Appendix.

Two important checks have been made in  to demonstrate that this is indeed the correct expression for the energy-momentum tensor. One is the consistency with the way the action of the noncommutative gauge theory  depends on the bulk metric $`g_{MN}`$.

(47)
$$\begin{array}{c}S=\frac{1}{g_{YM}^2}dx\sqrt{detG}[\frac{1}{4}G^{MN}G^{PQ}(F_{MP}\theta _{MP}^1)(F_{NQ}\theta _{NQ}^1)\hfill \\ \hfill +\frac{1}{2}G^{MN}g_{\alpha \beta }D_M\mathrm{\Phi }^\alpha D_N\mathrm{\Phi }^\beta \frac{1}{4}g_{\alpha \beta }g_{\gamma \delta }[\mathrm{\Phi }^\alpha ,\mathrm{\Phi }^\gamma ][\mathrm{\Phi }^\beta ,\mathrm{\Phi }^\delta ]].\end{array}$$

Here $`\theta _{MN}^1=0`$ unless $`(M,N)=(i,j)`$. Since the action is derived assuming that the bulk metric is flat, we can compare the variation $`S/g_{MN}`$ with the zero-momentum limit of the energy-momentum tensor.

The other one is the conservation of the energy-momentum tensor,

(48)
$$k_MT^{MN}=0.$$

We now have a simpler proof of this in the language of the matrix model discussed in the preceding section.

One may have expected that the energy-momentum tensor of the noncommutative theory should reduce to that of the commutative theory (theory defined on the commutative space) in the limit $`\theta ^{ij}0`$.

(49)
$$T_{comm}^{ij}=G^{ii^{}}G^{jj^{}}(F_{i^{}m}F_{j^{}n}G^{mn}\frac{1}{4}G_{i^{}j^{}}F^2+\mathrm{}),\text{etc}.$$

Here $`\mathrm{}`$ denotes terms containing the scalar field $`\mathrm{\Phi }^\alpha `$ $`etc`$. It turned out that they are different in the following three accounts:

1. In the commutative case (49), the metric $`G^{ii^{}}G^{jj^{}}`$ is used to raise the indices of the energy-momentum tensor, whereas in the noncommutative case (46) we use $`\theta ^{ii^{}}\theta ^{jj^{}}`$ in the corresponding term.
2. In the noncommutative case, the field strength $`F_{ij}`$ is shifted by $`\theta _{ij}^1`$.
3. In the noncommutative case, there are no terms that correspond to the term $`\frac{1}{4}G^{ij}F^2`$ in (49).

All these differences are important in order to maintain the consistency with the metric dependence of the Seiberg-Witten action (47) and the conservation of the energy-momentum tensor.

## 4. Bosonic string

We can use the same technique to derive the energy-momentum tensor of multiple D-branes in the bosonic string theory in the zero-slope limit. The closed string vertex operator is

(50)
$$V(z)=h_{MN}(k)X^M\overline{}X^Ne^{ikX(z)},$$

and the open string vertex operator is nothing but the bosonic part of $`U^{(0)}(t)`$ in (18). In the superstring case, the contribution to the graviton coupling starts at $`O(\alpha _{}^{}{}_{}{}^{3})`$. In the bosonic case, on the other hand, it starts at $`O(\alpha _{}^{}{}_{}{}^{2})`$ coming from the product of the following contractions:

(51) $`X^I(z)i_{}X^K(t)`$ $`={\displaystyle \frac{\alpha ^{}g^{IK}}{(zt)^2}}={\displaystyle \frac{2\pi i\alpha ^{}g^{IK}}{z\overline{z}}}e^{2\pi i\tau (t,z)}{\displaystyle \frac{\tau (t,z)}{t}},`$
$`\overline{}X^J(\overline{z})i_{}X^L(t)`$ $`={\displaystyle \frac{\alpha ^{}g^{JL}}{(\overline{z}t)^2}}={\displaystyle \frac{2\pi i\alpha ^{}g^{JL}}{z\overline{z}}}e^{2\pi i\tau (t,z)}{\displaystyle \frac{\tau (t,z)}{t}}.`$
The expression for $`T^{IJ}`$ is then

(52) $`T^{IJ}`$ $`={\displaystyle }dxe^{ikx}tr[\mathrm{exp}\left(i{\displaystyle _0^1}d\tau k_IY^I(x)\right)`$
$`\times [{\displaystyle _0^1}d\tau _1e^{2\pi i\tau _1}Y^I(x){\displaystyle _0^1}d\tau _2e^{2\pi i\tau _2}Y^J(x)+(IJ)]]`$
$`={\displaystyle }dxe^{ikx}tr{\displaystyle _0^1}d\tau [e^{2\pi i\tau }Y^Ie^{i\tau k_IY^I}Y^Je^{i(1\tau )k_IY^I}+(IJ)].`$
It is interesting to note that this leading term in the zero-slope limit vanishes in the case of a single brane since all the scalar fields commute and the only $`\tau `$ dependence of the integrand is $`e^{2\pi i\tau }`$. Thus $`_0^1𝑑\tau e^{2\pi i\tau }=0`$. Another consequence of the existence of the factor $`e^{2\pi i\tau }`$ is that the ordering of the scalar fields no longer obeys the symmetrized trace prescription,

(53)
$$T^{IJ}2𝑑xe^{ikx}Str\left[e^{ik_IY^I}Y^IY^J\right].$$

One can see this easily by considering the case of a single D-brane where $`T^{IJ}=0`$ in this limit as we just mentioned whereas the symmetrized trace on the right-hand side is non-vanishing.

To see the relation of this energy-momentum tensor to those of the noncommutative gauge theories derived in  (reproduced in the Appendix of this paper), the following identity is useful, which can be derived by integration by parts,

(54) $`tr{\displaystyle _0^1}𝑑\tau e^{\pm 2\pi i\tau }Ae^{i\tau kY}Be^{i(1\tau )kY}`$
$`=tr{\displaystyle _0^1}𝑑\tau {\displaystyle \frac{1}{\pm 2\pi i}}{\displaystyle \frac{d}{d\tau }}\left(e^{\pm 2\pi i\tau }\right)Ae^{i\tau kY}Be^{i(1\tau )kY}`$
$`={\displaystyle \frac{1}{\pm 2\pi i}}tr[A,B]e^{ikY}`$
$`{\displaystyle \frac{1}{\pm 2\pi i}}tr{\displaystyle _0^1}𝑑\tau e^{\pm 2\pi i\tau }Ae^{i\tau kY}[ikY,B]e^{i(1\tau )kY}`$
$`={\displaystyle \frac{1}{\pm 2\pi i}}tr[A,B]e^{ikY}`$
$`{\displaystyle \frac{1}{(2\pi )^2}}tr[A,[ikY,B]]e^{ikY}`$
$`+{\displaystyle \frac{1}{(2\pi )^2}}tr{\displaystyle _0^1}𝑑\tau e^{\pm 2\pi i\tau }[ikY,A]e^{i\tau kY}[ikY,B]e^{i(1\tau )kY}.`$
After we symmetrize with respect to $`A`$ and $`B`$, the first term on the right-hand side vanishes. Combining this with

$`[ik_IY^I,x^j+\theta ^{jm}A_m]=\theta ^{jm}\left(l^n(F_{mn}\theta _{mn}^1)+y_\beta D_m\mathrm{\Phi }^\beta \right),`$
$`[x^i+\theta ^{im}A_m,𝒪]=i\theta ^{im}D_m𝒪,`$
it is then straightforward to see that (52) reproduces the energy-momentum tensors of the noncommutative theories (62)–(64).

The proof of the conservation of the energy-momentum tensor $`k_IT^{IJ}=0`$ in the matrix form (52) of $`T^{IJ}`$ is given by

(55) $`tr{\displaystyle _0^1}𝑑\tau e^{\pm 2\pi i\tau }kYe^{i\tau kY}Y^Je^{i(1\tau )kY}`$
$`=tr{\displaystyle _0^1}𝑑\tau e^{\pm 2\pi i\tau }e^{i\tau kY}kYY^Je^{i(1\tau )kY}`$
$`={\displaystyle _0^1}𝑑\tau e^{\pm 2\pi i\tau }trkYY^Je^{ikY}`$
$`=0\text{because of the }\tau \text{ integration}.`$
Note that the conservation holds without using the equation of motion. This proof is significantly simpler than the one given in .

## Acknowledgments

We would like to thank Sumit Das, Satoshi Iso, Washington Taylor and Sandip Trivedi for discussions. H. O. thanks the organizers of the conference, Strings 2001, at Tata Institute of Fundamental Research, Mumbai, India for their hospitality and for giving him the opportunity to present our work at the conference.

The research was supported in part by the DOE grant DE-FG03-92ER40701 and the Caltech Discovery Fund.

## Appendix A Energy-momentum tensors
of noncommutative gauge theories

### A.1. Superstring

(56)
$$\begin{array}{c}T^{ij}=2\alpha ^{}\theta ^{ii^{}}\theta ^{jj^{}}\hfill \\ \hfill \times dx[e^{ikx}\mathrm{exp}\left(i_0^1d\tau (l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\right)\\ \hfill \times _0^1d\tau _1_0^1d\tau _2[G^{mn}(F_{i^{}m}(x+l\tau _1)\theta _{i^{}m}^1)(F_{j^{}n}(x+l\tau _2)\theta _{j^{}n}^1)\\ \hfill +G^{\mu \nu }F_{i^{}\mu }(x+l\tau _1)F_{j^{}\nu }(x+l\tau _2)+g_{\alpha \beta }D_i^{}\mathrm{\Phi }^\alpha (x+l\tau _1)D_j^{}\mathrm{\Phi }^\beta (x+l\tau _2)]],\end{array}$$

(57)
$$\begin{array}{c}T^{\mu j}=2\alpha ^{}G^{\mu \mu ^{}}\theta ^{jj^{}}\hfill \\ \hfill \times dx[e^{ikx}\mathrm{exp}\left(i_0^1d\tau (l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\right)\\ \hfill \times _0^1d\tau ^{}F_{\mu ^{}j^{}}(x+l\tau ^{})],\end{array}$$

(58)
$$\begin{array}{c}T^{\alpha j}=4\pi \alpha _{}^{}{}_{}{}^{2}\theta ^{jj^{}}\hfill \\ \hfill \times dx[e^{ikx}\mathrm{exp}\left(i_0^1d\tau (l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\right)\\ \hfill \times _0^1d\tau _1_0^1d\tau _2[G^{mn}D_m\mathrm{\Phi }^\alpha (x+l\tau _1)(F_{j^{}n}(x+l\tau _2)\theta _{j^{}n}^1)\\ \hfill +G^{\mu \nu }D_\mu \mathrm{\Phi }^\alpha (x+l\tau _1)F_{j^{}\nu }(x+l\tau _2)\\ \hfill +g_{\beta \gamma }i[\mathrm{\Phi }^\beta ,\mathrm{\Phi }^\alpha ](x+l\tau _1)D_j^{}\mathrm{\Phi }^\gamma (x+l\tau _2)]],\end{array}$$

(59)
$$\begin{array}{c}T^{\mu \nu }=2\alpha ^{}G^{\mu \nu }\hfill \\ \hfill \times dx\left[e^{ikx}\mathrm{exp}\left(i_0^1d\tau (l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\right)\right],\end{array}$$

(60)
$$\begin{array}{c}T^{\alpha \nu }=4\pi \alpha _{}^{}{}_{}{}^{2}G^{\nu \nu ^{}}\hfill \\ \hfill \times dx[e^{ikx}\mathrm{exp}\left(i_0^1d\tau (l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\right)\\ \hfill \times _0^1d\tau ^{}D_\nu ^{}\mathrm{\Phi }^\alpha (x+l\tau ^{})],\end{array}$$

(61)
$$\begin{array}{c}T^{\alpha \beta }=8\pi ^2\alpha _{}^{}{}_{}{}^{3}\hfill \\ \hfill \times dx[e^{ikx}\mathrm{exp}\left(i_0^1d\tau (l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\right)\\ \hfill \times _0^1d\tau _1_0^1d\tau _2[G^{ij}D_i\mathrm{\Phi }^\alpha (x+l\tau _1)D_j\mathrm{\Phi }^\beta (x+l\tau _2)\\ \hfill +G^{\mu \nu }D_\mu \mathrm{\Phi }^\alpha (x+l\tau _1)D_\nu \mathrm{\Phi }^\beta (x+l\tau _2)\\ \hfill g_{\gamma \delta }[\mathrm{\Phi }^\gamma ,\mathrm{\Phi }^\alpha ](x+l\tau _1)[\mathrm{\Phi }^\delta ,\mathrm{\Phi }^\beta ](x+l\tau _2)]].\end{array}$$

### A.2. Bosonic string

(62)
$$\begin{array}{c}T^{ij}=\frac{\theta ^{ii^{}}\theta ^{jj^{}}+\theta ^{ji^{}}\theta ^{ij^{}}}{(2\pi )^2}\hfill \\ \hfill \times dx[e^{ikx}\mathrm{exp}(i_0^1d\tau l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\\ \hfill \times \{i_0^1d\tau _1e^{2\pi i\tau _1}(l^mF_{i^{}m}(x+l\tau _1)+y_\alpha D_i^{}\mathrm{\Phi }^\alpha (x+l\tau _1))\\ \hfill \times i_0^1𝑑\tau _2e^{2\pi i\tau _2}\left(l^nF_{j^{}n}(x+l\tau _2)+y_\beta D_j^{}\mathrm{\Phi }^\beta (x+l\tau _2)\right)\\ \hfill +i_0^1d\tau ^{}(l^mD_i^{}F_{j^{}m}(x+l\tau ^{})+y_\alpha D_i^{}D_j^{}\mathrm{\Phi }^\alpha (x+l\tau ^{}))\}].\end{array}$$

(63)
$$\begin{array}{c}T^{\alpha j}=\alpha ^{}\theta ^{jm}dx[e^{ikx}\mathrm{exp}(i_0^1d\tau l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\hfill \\ \hfill \times \{_0^1d\tau _1e^{2\pi i\tau _1}\mathrm{\Phi }^\alpha (x+l\tau _1)\\ \hfill \times i_0^1𝑑\tau _2e^{2\pi i\tau _2}\left(l^nF_{mn}(x+l\tau _2)+y_\beta D_m\mathrm{\Phi }^\beta (x+l\tau _2)\right)\\ \hfill _0^1𝑑\tau _1e^{2\pi i\tau _1}\mathrm{\Phi }^\alpha (x+l\tau _1)\\ \hfill \times i_0^1d\tau _2e^{2\pi i\tau _2}(l^nF_{mn}(x+l\tau _2)+y_\beta D_m\mathrm{\Phi }^\beta (x+l\tau _2))\}].\end{array}$$

(64)
$$\begin{array}{c}T^{\alpha \beta }=(2\pi \alpha ^{})^2\hfill \\ \hfill \times dx[e^{ikx}\mathrm{exp}(i_0^1d\tau l^iA_i(x+l\tau )+y_\alpha \mathrm{\Phi }^\alpha (x+l\tau ))\\ \hfill \times \{_0^1d\tau _1e^{2\pi i\tau _1}\mathrm{\Phi }^\alpha (x+l\tau _1)_0^1d\tau _2e^{2\pi i\tau _2}\mathrm{\Phi }^\beta (x+l\tau _2)+(\alpha \beta )\}].\end{array}$$


format: text

text:
# The Nystrom plus Correction Method for Solving Bound State Equations in Momentum Space

## I Introduction

The momentum space Schrödinger equation has a singular kernel for both the Coulomb and linear potentials. The Coulomb singularity is removed with the Lande subtraction method . Previous work  showed how to remove the singularity from the linear potential using a subtraction method with basis functions. A problem with this method is that one must guess a suitable set of basis functions in advance. In this paper, we show that Nystrom method  can solve the same problem more simply amd accurately. We begin with a review of the basis function method. Then we introduce the Nystrom method and apply it to the s-state momentum space Schrödinger equation with a linear potential. We use our new numerical results to show that the Nystrom plus correction method is more accurate than the basis function method. At the end, we generalize the Nystrom method to higher angular momentum quantum numbers ($`l>0`$).

## II The Basis Function Method

We begin this paper with a discussion of the basis function method to give the proper theoretical motivation. We shall use the simplest momentum space Schrödinger equation to illustrate the principles of the numerical methods, which is the s-state equation.

The momentum space Schrödinger equation is related to an integral equation of the Fredholm type

$$K(p,p^{})\varphi (p^{})𝑑p^{}=\lambda \varphi (p).$$
(1)
Suppose that the wavefunction $`\varphi `$ can be expanded in a set of basis functions $`\{g_i\}`$, such that

$$\varphi (p)=\underset{i=1}{\overset{N}{}}C_ig_i(p),$$
(2)
where $`C_i`$ are constant coefficients. Substitute Eq.  into Eq.  to obtain

$$\underset{i=1}{\overset{N}{}}K(p,p^{})C_ig_i(p^{})𝑑p^{}=\lambda \underset{i=1}{\overset{N}{}}C_ig_i(p).$$
(3)
Now multiply both sides of Eq.  with $`g_j(p)`$ and integrate over $`p`$ to symmetrize the equation over $`i`$ and $`j`$,

$$\underset{i=1}{\overset{N}{}}C_i\underset{A_{ij}}{\underset{}{K(p,p^{})g_i(p^{})g_j(p)𝑑p^{}𝑑p}}=\lambda \underset{i=1}{\overset{N}{}}C_i\underset{B_{ij}}{\underset{}{g_i(p)g_j(p)𝑑p}},$$
(4)
and the result is a matrix equation,

$$\underset{i=1}{\overset{N}{}}A_{ij}C_i=\lambda \underset{i=1}{\overset{N}{}}B_{ij}C_i,$$
(5)
where $`C_i`$ is the eigenvector and $`\lambda `$ the eigenvalue. The indices $`i`$ and $`j`$ correspond to the quadrature points $`p`$ and $`p^{}`$. $`N`$ represents the number of mesh points. In the case of the momentum space Schrödinger equation with a Coulomb or linear potential, the kernel $`𝐀`$ is singular. A simple example is the momentum space Schrödinger equation with a linear potential in the $`S`$-state ,

$$\frac{p^2}{2\mu }\varphi _{n0}(p)+\frac{\lambda _L}{\pi p^2}_0^{\mathrm{}}\underset{V_0^L(p,p^{})}{\underset{}{\left[\frac{\eta ^2}{p^{}p}Q_0^{\prime \prime }(y)+Q_0^{}(y)\right]}}\varphi _{n0}(p^{})𝑑p^{}=E_{n0}\varphi _{n0}(p),$$
(6)
where $`y=(p^2+p^2)/2p^{}p`$,

$$Q_0^{}(y)=p^{}p\left[\frac{1}{(p+p^{})^2+\eta ^2}\frac{1}{(pp^{})^2+\eta ^2}\right],$$
(7)
and

$$\frac{\eta ^2}{p^{}p}Q_0^{\prime \prime }(y)=\eta ^2\left(p^2+p^2+\eta ^2\right)\left[\frac{1}{(p+p^{})^2+\eta ^2}\frac{1}{(pp^{})^2+\eta ^2}\right]^2.$$
(8)
Lande subtraction  involves subtracting a zero term

$$_0^{\mathrm{}}\left[\frac{\eta ^2}{p^{}p}Q_0^{\prime \prime }(y)+Q_0^{}(y)\right]𝑑p^{}=0$$
(9)
from Eq.  such that

$$\frac{p^2}{2\mu }\varphi _{n0}(p)+\frac{\lambda _L}{\pi p^2}_0^{\mathrm{}}\left[\frac{\eta ^2}{p^{}p}Q_0^{\prime \prime }(y)+Q_0^{}(y)\right][\varphi _{n0}(p^{})\varphi _{n0}(p)]𝑑p^{}=E_{n0}\varphi _{n0}(p).$$
(10)
Using Eqs. , the integral in Eq.  for $`p>0`$ in the limit of $`y1`$ can be shown to equal

$$\underset{\eta 0}{lim}\underset{pp^{}}{lim}\frac{\lambda _L}{\pi }\left[2\eta ^2\left(\frac{1}{(pp^{})^2+\eta ^2}\right)^2\frac{1}{(pp^{})^2+\eta ^2}\right](pp^{})^2\frac{d\varphi _{n0}}{dp}=0.$$
(11)
The order of the limits in Eq.  is important. The reverse order will lead to the nonsensical result $`Q_0^{}(y)𝑑p^{}=0`$. Next, in the limit of $`p,p^{}0`$, $`(p+p^{})^2=(pp)^2`$. By substituting this equality into Eqs. , it can be shown again that the integral in Eq.  vanishes for $`p0`$ at $`y=1`$. At the end, the integral vanishes at $`y=1,p`$. Away from the singularities, both integrands in the integral of Eq.  are finite. By taking $`\eta 0`$, the first integrand vanishes. The final form of Eq.  is

$$\frac{p^2}{2\mu }\varphi _{n0}(p)+\frac{\lambda _L}{\pi p^2}_0^{\mathrm{}}Q_0^{}(y)[\varphi _{n0}(p^{})\varphi _{n0}(p)]𝑑p^{}=E_{n0}\varphi _{n0}(p),$$
(12)
where $`Q_0^{}(y)=1/(1y^2)`$. As mentioned before, $`\varphi `$ is expanded in basis functions, followed by integrating Eq.  over $`p`$ to generate a matrix equation. The basis functions used in previous publications  are

$$g_i^A(p)=exp\left[\frac{p^2i^2}{M}\right]$$
(13)
and

$$g_i^B(p)=\frac{1}{(i/M)^2+p^4}$$
(14)
respectively, where $`M`$ is the maximum number of basis functions used. $`M`$ has a maximum because the code crashes when too many basis functions are used. The basis functions $`g_i^A(p)`$ and $`g_i^B(p)`$ have $`M=18`$.

The singularity of the kernel is a major challenge in solving the integral equation with a linear potential. It was shown  that a simple pole remains even after subtraction. The role of the basis functions is to make possible the evaluation of the Cauchy principal value of the subtracted integral using the Sloan method . To illustrate the Sloan method, we suppose that $`f(x)`$ has a simple pole such that

$$f(x)=\frac{g(x)}{x}$$
(15)
where $`g(x)`$ is regular. The Cauchy principal value of the subtracted integral of $`f(x)`$ can be evaluated if the range of integration is symmetric. For example, the numerical integration of

$$_1^1\frac{g(x)g(0)}{x}𝑑x$$
(16)
yields the Cauchy principal value because the point $`x=0`$ is skipped when quadrature points are generated in the symmetric interval $`(1,1)`$. The subtraction term has zero contribution since

$$_1^1\frac{1}{x}𝑑x=0.$$
(17)
The purpose of this term is to justify the existence of the Cauchy principle value and to reduce numerical errors. In order to apply the Sloan method to Eq. , the integration variable is transformed from $`p`$ to $`x`$ such that $`x`$ is centered at zero and its range is symmetric.

In the case of the Coulomb potential, the kernel has a logarithmic singularity,

$$Q_0(y)=\frac{1}{2}\mathrm{ln}\left|\frac{y+1}{y1}\right|,$$
(18)
which is completely removed by Lande subtraction  because no simple pole remains after the subtraction.

The key to the success of the basis function method is the availability of a suitable set of basis functions for a particular problem. Unfortunately there is no *a priori* reason why the same set of basis functions will work in every situation. For this reason, it may be advantageous to have a method (such as the Nystrom method) that does not depend on a choice of the basis functions.

## III The Nystrom Method

In general, an integral equation of the Fredholm type

$$G(p)\varphi (p)+_0^{\mathrm{}}F(p,p^{})\varphi (p^{})𝑑p^{}=\lambda \varphi (p)$$
(19)
can be rewritten as a matrix equation as

$$\underset{j=1}{\overset{N}{}}K_{ij}\varphi _j\underset{j=1}{\overset{N}{}}\left(G_{ii}+F_{ij}\right)\varphi _j=\lambda \varphi _i,$$
(20)
where $`K_{ij}`$ is the kernel and $`i`$ and $`j`$ are now indices corresponding to $`p`$ and $`p^{}`$. Instead of integrating over $`p`$ from 0 to $`\mathrm{}`$, we integrate over $`x`$ from $`1`$ to 1. Transform $`x_i`$ to $`p_i`$ by the transformation

$$p(x)=\mathrm{tan}\left(\frac{1+x}{4}\pi \right).$$
(21)
The mesh points $`x_i`$ and the weights $`wt_i`$ are generated by the gaussian quadrature rule using the routine gauleg from Numerical Recipes . In order to integrate along $`x[1,1]`$ instead of $`p[0,\mathrm{})`$, Eq.  is transformed as

$$G(x)\varphi (x)+_1^1F(x,x^{})\varphi (x^{})\frac{dp^{}}{dx^{}}𝑑x^{}=\lambda \varphi (x)$$
(22)
Changing the dummy variable inside the integral and substituting the differentiation of Eq.  with

$$dp=\frac{\pi }{4}\mathrm{sec}^2\left(\frac{1+x}{4}\pi \right)dx=\frac{\pi }{4}\left(1+p^2\right)dx$$
(23)
gives

$$\frac{p_i^2}{2\mu }\varphi _i+\frac{\lambda _L}{4p_i^2}_1^1Q_0^{}(y1)[\varphi _j\varphi _i]\mathrm{sec}^2\left(\frac{1+x_j}{4}\pi \right)𝑑x_j=E\varphi _i.$$
(24)
Eq.  can now be written as a matrix equation,

$`{\displaystyle \frac{p_i^2}{2\mu }}\varphi _i+{\displaystyle \frac{\lambda _L}{4p_i^2}}{\displaystyle \underset{j=1}{\overset{N}{}}}Q_0^{}(y1)\varphi _j\mathrm{sec}^2\left({\displaystyle \frac{1+x_j}{4}}\pi \right)wt_j`$ (25)
$`{\displaystyle \frac{\lambda _L}{4p_i^2}}\varphi _i{\displaystyle \underset{k=1}{\overset{N}{}}}Q_0^{}(y1)\mathrm{sec}^2\left({\displaystyle \frac{1+x_k}{4}}\pi \right)wt_k=E\varphi _i.`$ (26)
The left hand side of Eq.  is the kernel times the eigenvector and the right hand side is the eigenvalue times the eigenvector. The sum over $`k`$ is independent of the eigenvector, which is just a scalar. The terms on the left that have only one index $`i`$ belong to the diagonal elements $`K_{ii}`$. The terms with mixed indices $`i`$ and $`j`$ make up the off-diagonal elements, $`K_{ij}`$. More explicitly, the matrix elements of the kernel are

$`K_{ii}`$ $`=`$ $`{\displaystyle \frac{p_i^2}{2\mu }}{\displaystyle \frac{\lambda _L}{4p_i^2}}{\displaystyle \underset{k}{}}Q_0^{}(y1)\left(1+p_k^2\right)wt_k,`$ (27)
$`K_{ij}`$ $`=`$ $`{\displaystyle \frac{\lambda _L}{4p_i^2}}Q_0^{}(y)\left(1+p_j^2\right)wt_j,(ij).`$ (28)
So far the kernel $`𝐊`$ is asymmetric under the interchange of $`i`$ and $`j`$. We can improve the stability and the efficiency of the numerical solutions by symmetrizing Eq. . We do so by multiplying the equation with $`p_i^2\left(1+p_i^2\right)`$. It will change the original matrix equation

$$𝐊𝐱=\lambda 𝐱$$
(29)
to an equivalent matrix equation

$$𝐊^{}𝐱=\lambda 𝐂𝐱,$$
(30)
where $`𝐂`$ is a diagonal matrix and $`𝐊^{}=𝐂𝐊`$. If $`𝐂`$ is positive definite, meaning

$$𝐱^T𝐂𝐱0,\mathrm{vectors}𝐱,$$
(31)
then $`𝐂`$ can be Cholesky-decomposed as

$$𝐂=𝐋𝐋^T,$$
(32)
where $`𝐋`$ is a unique lower triangular matrix. The reason for performing Cholesky decomposition is that the new matrix

$$𝐊^{\prime \prime }𝐋^1𝐊^{}(𝐋^1)^T$$
(33)
is real symmetric and yields the same eigenvalues as Eq. . In the case of $`C_{ii}=p_i^2\left(1+p_i^2\right)`$, $`𝐂`$ is guaranteed to be positive definite. After symmetrization, Eq.  does not change ($`K_{ii}^{\prime \prime }=K_{ii}`$) while Eq.  becomes

$`K_{ij}^{\prime \prime }`$ $`=`$ $`{\displaystyle \frac{\lambda _L}{4p_ip_j}}Q_0^{}(y)\sqrt{\left(1+p_i^2\right)\left(1+p_j^2\right)wt_iwt_j},(ij).`$ (34)
The eigenvalues of $`𝐊^{\prime \prime }`$ can be calculated by using standard packages such as EISPAK. In this paper, we use the tred2 and tqli routines in *Numerical Recipes* .

## IV The Correction Method

Maung, Norbury and Kahana  have shown that the subtraction method does not completely remove the singularity at $`y=1`$. There is a residual simple pole term

$$\frac{4\lambda _L}{\pi }\frac{d\varphi _{n0}}{dp}_0^{\mathrm{}}\frac{p^2}{(p^{}+p)^2(p^{}p)}𝑑p^{}$$
(35)
remaining after the subtraction. The basis function method evaluates the Cauchy principal value by the Sloan method as described in Section 2. The Sloan method eliminates the simple pole term by integrating symmetrically around the singularity. Symmetrical integration involves splitting the range of integration is into two intervals,

$$_0^{\mathrm{}}𝑑p^{}=_0^{2p}𝑑p^{}+_{2p}^{\mathrm{}}𝑑p^{}.$$
(36)
The singularity at $`p=p^{}`$ is contained in the first term on the right hand side of Eq.  which is assigned a symmetric transformation rule $`(dp^{}/dx)_1`$. The second term generally has a different transformation rule $`(dp^{}/dx)_2`$ because it is mapping between two different sets, namely $`(2p,\mathrm{})(1,M]`$ (for some real number $`M`$), such that

$$_0^{\mathrm{}}𝑑p^{}_1^1\left(\frac{dp^{}}{dx^{}}\right)_1𝑑x^{}+_1^M\left(\frac{dp^{}}{dx^{}}\right)_2𝑑x^{}.$$
(37)
Notice that the division of the range of integration moves with $`p`$. If two transformation rules are used with a moving division, each row (column) of the kernel has a different way to map $`[0,\mathrm{})`$ to $`[1,M]`$. But the eigenvector $`\varphi (p)`$ must be mapped to $`\varphi _i`$ in a unique way. This mismatch between the mappings of the kernel and the eigenvector does not affect the basis function method (see Eq. \[2.24\] of reference )

$`{\displaystyle \underset{i=1}{\overset{N}{}}}C_i\left[{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{p^4}{2\mu }}g_j(p)g_i𝑑p+{\displaystyle \frac{\lambda _L}{\pi }}{\displaystyle _0^{\mathrm{}}}{\displaystyle _0^{\mathrm{}}}Q_0^{}(y)g_j(p)\left[g_i(p^{})g_i(p)\right]𝑑p^{}𝑑p\right]`$ (38)
$`=`$ $`E{\displaystyle \underset{i=1}{\overset{N}{}}}C_i{\displaystyle _0^{\mathrm{}}}p^2g_j(p)g_i(p)𝑑p`$ (39)
because the eigenvector $`C_i`$ is an $`N`$-tuple of the coefficients of the basis function expansion of the wavefunction $`\varphi (p)`$ and is independent of the transformation rules. In the case of the Nystrom method, the problem is real, at least for the range of integration that we are interested in. Therefore we cannot evaluate the Cauchy principal value by symmetric integration in the Nystrom method. In other words, a new method is needed to treat the errors arising from the simple pole term.

So far the error term Eq.  is not included in the Nystrom kernel in our derivation and is contributing to the errors of the eigenvalues. Since the error term Eq.  involves $`d\varphi /dp`$, we associate it with the error of the wavefunction

$$\mathrm{\Delta }\varphi =\mathrm{\Delta }p\frac{d\varphi }{dp}\frac{1}{N}\frac{d\varphi }{dp},$$
(40)
where the mesh size $`\mathrm{\Delta }p`$ has an $`N^1`$ dependence. This fact leads to an estimate of the $`N`$ dependence of the error of the eigenvalue, $`\mathrm{\Delta }E`$. Let the approximate eigenvector to be $`\varphi ^{}`$ and the approximate eigenvalue $`E^{}`$. It is reasonable to say that an approximate kernel $`K`$ acting on an approximate eigenvector $`\varphi ^{}`$ yields an approximate eigenvalue $`E^{}`$ as in

$`K\varphi ^{}=E^{}\varphi ^{}`$ (41)
$``$ $`K(\varphi +\mathrm{\Delta }\varphi )=(E+\mathrm{\Delta }E)(\varphi +\mathrm{\Delta }\varphi ).`$ (42)
It is easy to see that

$`\mathrm{\Delta }E`$ $``$ $`(KE){\displaystyle \frac{\mathrm{\Delta }\varphi }{\varphi }}`$ (43)
$`=`$ $`\left({\displaystyle \frac{KE}{\varphi }}\right){\displaystyle \frac{d\varphi }{dp}}\mathrm{\Delta }p`$ (44)
$`=`$ $`ϵ{\displaystyle \frac{1}{N}}`$ (45)
It is safe to assume that $`(KE)d\varphi /dp<<1`$. $`\varphi ^1`$ can be interpreted as the normalization. The product of all of the pseudo-constants is labelled as the coefficient $`ϵ`$. The approximate eigenvalue $`E^{}`$ produced in the background of Eq.  is related to the true eigenvalue $`E`$ by

$$E_n^{}=E_n+ϵf_n(N),$$
(46)
where $`n`$ is the principal quantum number, $`ϵ`$ a constant and $`f_n(N)`$ is a function approximately equals to $`N^1`$. In general, $`f_n(N)`$ varies slightly depending on the type of integral equation and the potential. As a first order approximation, assume that

$$f_n(N)=N^{1\alpha (n1)}.$$
(47)
The exponent of Eq.  is a first order Taylor series expansion of some negative unity function around $`n=1`$. The constant $`\alpha `$ is always taken to be small. More particularly, choose an $`\alpha `$ such that the variance of $`E_n`$ and $`ϵ`$ and $`\chi ^2`$ are minimized in the linear fit. Finally the refinement of an eigenvalue involves generating a set of $`E_n^{}`$ for various $`N`$ by the Nystrom method and then extrapolating $`E_n`$ by a $`\chi ^2`$ linear fit in the graph of $`E_n`$ versus $`f_n(N)`$. In the case of Eq. , $`\alpha =0.004`$ is an optimal choice. The numerical results are explained in Section 6.

The order of the Nystrom algorithm is derived from those of tred2 and tqli, which is $`𝒪(N^2)`$ , compared with the basis function’s $`𝒪(M^2N)`$, which comes from the product of the size of the matrix $`M^2`$ and the number of integration mesh points. $`N`$ is typically around 1000 and $`M`$ is 20. The basis function method is generally more efficient than the Nystrom method. However, for any given set of basis functions, the accuracy of the eigenvalues cannot be improved arbitrarily by increasing the number of basis functions because $`M`$ is bounded from above due to numerical errors. The prospect of improving the accuracy of the basis function algorithm depends on the availability of a set of more suitable basis functions for a specific problem. In the case of the Nystrom plus correction method, accuracy is optimized automatically by the correction scheme. The numerical results obtained by the Nystrom and basis function methods are quoted with optimal accuracy in this paper.

## V Exact S-state Solution

The eigenvalue of Eq.  can be solved exactly in configuation space. We shall use the analytic results to check our numerical results. The non-relativistic Schrödinger equation can be written as

$$\left(\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}\right)R2\mu [\lambda _LrE]R=0.$$
(48)
Let $`SrR`$, then Eq.  can be simplified as

$$\frac{d^2}{dr^2}S2\mu [\lambda _LrE]S=0.$$
(49)
Define a new variable

$$x\left(\frac{2\mu }{\lambda _L^2}\right)^{\frac{1}{3}}[\lambda _LrE],$$
(50)
such that Eq.  can be transformed as

$$S^{\prime \prime }xS=0,$$
(51)
which is the Airy equation. The solution which satisfies the boundary condition $`S0`$ as $`x\mathrm{}`$ is the Airy function $`\mathrm{Ai}(x)`$. It is easy to show that the eigen-energy formula is

$$E_n=x_n\left(\frac{\lambda _L^2}{2\mu }\right)^{\frac{1}{3}},$$
(52)
where $`x_n`$ is the $`n`$-th zero of the Airy function counting from $`x=0`$ along $`x`$. In reference , the values $`\lambda _L=5`$ and $`\mu =0.75`$ are used. In this case, the eigen-energy formula is

$$E_n=2.554364772x_n.$$
(53)
## VI Numerical Results for the S-state

The accuracy of the Nystrom plus correction method is sensitive to the range of $`N`$. In this paper, increments of 100 in the range of $`100N1400`$ are used. The reason for this choice is that there are not enough spacings between the eigenvalues for $`N<100`$ and for $`N>1500`$ the numerical noise begins to corrupt the monotonic convergent behavior of the eigenvalues. The correct eigenvalues are extrapolated from these numerical data by a $`\chi ^2`$ linear fit as described in Section 4. The exact $`S`$-state eigenvalues are tabulated against the numeical results obtained by the basis function method and the Nystrom plus corrections method in Table \[I\]. It shows that the numerical results obtained by the Nystrom method plus corrections are more accurate than the results obtained by the basis function method.

The kernel written for the non-relativistic Schrödinger equation can be easily generalized to that of the relativistic 2-body Thompson equation in the center-of-mass frame by the replacement

$$\frac{p^2}{2\mu }2\left(\sqrt{p^2+m^2}m\right),$$
(54)
where $`\mu `$ is the reduced mass and $`m`$ is the mass of each of the two equal mass elementary particles. The numerical results obtained using the Thompson equation is compared against those using the non-relativistic Schrödinger equation in Table \[II\] calculating to 2 decimal places. Our new results are exactly the same as the previous results obtained in reference  that uses basis functions $`g_i^A(p)`$ from Eq. .

## VII $`l0`$ Kernels

The $`l0`$ kernels for the linear and Coulomb potentials contain the Legendre function of the second kind $`Q_l(y)`$ and its derivative respectively. There are several mathematical issues that need to be addressed before constructing the $`l0`$ kernels. First of all, the definition of

$$y\frac{p^2+p^2}{2p^{}p}=\frac{1}{2}\left(\frac{p}{p^{}}+\frac{p^{}}{p}\right)$$
(55)
is easily seen to yield $`y>1`$ for $`p,p^{}>0`$. In reference , Maung *et al.* use the Legendre identity

$`Q_l(y)`$ $`=`$ $`P_l(y)Q_0(y)w_{l1}(y),`$ (56)
$`w_{l1}(y)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{l}{}}}{\displaystyle \frac{1}{m}}P_{lm}(y)P_{m1}(y),`$ (57)
which is valid for $`1y1`$  but can be extended to $`y>1`$ by analytic continuation . $`Q_l^{}(y)`$ is easily obtained by straightforward differentiation. The derivative of Legendre polynomial can be calculated from one of the recurrence formulas,

$$\frac{dP_l(y)}{dy}=y\frac{dP_{l1}(y)}{dy}+lP_{l1}(y),$$
(58)
which can be computed numerically by a recursive call. The Legendre function can be generated by modifying the routine plgndr in Numerical Recipes  to allow $`y>1`$. The accuracy of Eq.  and its derivatives are generally sufficient. Slightly more accurate results can be obtained by the explicit evaluation of the Neumann integral,

$$Q_l(y)=\frac{1}{2}_1^1\frac{1}{(yt)}P_l(t)𝑑t,$$
(59)
with derivative

$$Q_l^{}(y)=\frac{1}{2}_1^1\frac{1}{(yt)^2}P_l(t)𝑑t.$$
(60)
The first few $`Q_l(y)`$ are

$`Q_0(y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{y+1}{y1}},`$ (61)
$`Q_1(y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}y\mathrm{ln}{\displaystyle \frac{y+1}{y1}}1,`$ (62)
$`Q_2(y)`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(3y^21\right)\mathrm{ln}{\displaystyle \frac{y+1}{y1}}{\displaystyle \frac{3}{2}}y,`$ (63)
$`Q_3(y)`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(5y^33y\right)\mathrm{ln}{\displaystyle \frac{y+1}{y1}}{\displaystyle \frac{5}{2}}y^2+{\displaystyle \frac{2}{3}},`$ (64)
$`Q_4(y)`$ $`=`$ $`{\displaystyle \frac{1}{16}}\left(35y^430y^2+3\right)\mathrm{ln}{\displaystyle \frac{y+1}{y1}}{\displaystyle \frac{35}{8}}y^3+{\displaystyle \frac{55}{24}}y,`$ (65)
$`Q_5(y)`$ $`=`$ $`{\displaystyle \frac{1}{16}}\left(63y^570y^3+15y\right)\mathrm{ln}{\displaystyle \frac{y+1}{y1}}{\displaystyle \frac{63}{8}}y^4+{\displaystyle \frac{49}{8}}y^2{\displaystyle \frac{8}{15}}.`$ (66)
$`Q_l^{}(y)`$ can be obtained by the direct differentiation of $`Q_l(y)`$, such that

$`Q_0^{}(y)`$ $`=`$ $`{\displaystyle \frac{1}{1y^2}},`$ (67)
$`Q_1^{}(y)`$ $`=`$ $`{\displaystyle \frac{y}{1y^2}}{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{y1}{y+1}},`$ (68)
$`Q_2^{}(y)`$ $`=`$ $`{\displaystyle \frac{1}{1y^2}}{\displaystyle \frac{3}{2}}y\mathrm{ln}{\displaystyle \frac{y1}{y+1}}3,`$ (69)
$`Q_3^{}(y)`$ $`=`$ $`{\displaystyle \frac{y}{1y^2}}{\displaystyle \frac{15y^23}{4}}\mathrm{ln}{\displaystyle \frac{y1}{y+1}}{\displaystyle \frac{15}{2}}y,`$ (70)
$`Q_4^{}(y)`$ $`=`$ $`{\displaystyle \frac{1}{1y^2}}{\displaystyle \frac{35y^315y}{4}}\mathrm{ln}{\displaystyle \frac{y1}{y+1}}{\displaystyle \frac{35}{2}}y^2+{\displaystyle \frac{5}{3}},`$ (71)
$`Q_5^{}(y)`$ $`=`$ $`{\displaystyle \frac{y}{1y^2}}{\displaystyle \frac{315y^4210y^2+15}{16}}\mathrm{ln}{\displaystyle \frac{y1}{y+1}}{\displaystyle \frac{315}{8}}y^3+{\displaystyle \frac{105}{8}}y.`$ (72)
As $`y\mathrm{}`$, it is easily seen that $`Q_0(y)=Q_0^{}(y)0`$. This limit is true for all $`Q_l(y)`$ and $`Q_l^{}(y)`$ from applying the L’Hopital rule. Unfortunately straightforward numerical calculation of $`Q_l(y)`$ and $`Q_l^{}(y)`$ by using Eqs. \[61-72\] leads to serious numerical errors as $`y\mathrm{}`$. At the same time, it is observed that the numerical integration of Eq.  are reasonably accurate in the same regime. Therefore the two representations are combined to minimize numerical error by using the Neumann integrals for $`y>y_0`$ and the explicit formulas for $`yy_0`$. Our codes use $`y_0=50`$.

The subtracted momentum space NRSE with a linear potential is given in reference , which can be simplified as

$`{\displaystyle \frac{p^2}{2\mu }}\varphi _{nl}(p)+{\displaystyle \frac{\lambda _L}{\pi p^2}}{\displaystyle _0^{\mathrm{}}}Q_l^{}(y)\varphi _{nl}(p^{})𝑑p^{}{\displaystyle \frac{\lambda _L}{\pi p^2}}\varphi _{nl}(p){\displaystyle _0^{\mathrm{}}}Q_0^{}(y)𝑑p^{}`$ (73)
$`{\displaystyle \frac{\lambda _L}{\pi p}}{\displaystyle \frac{l(l+1)}{2}}\varphi _{nl}(p){\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{Q_0}{p^{}}}𝑑p^{}+{\displaystyle \frac{\lambda _L\pi }{p}}{\displaystyle \frac{l(l+1)}{4}}\varphi _{nl}(p)=E_{nl}\varphi _{nl}(p).`$ (74)
The matrix elements of a symmetric kernel for arbitrary $`l`$ are

$`K_{ii}`$ $`=`$ $`{\displaystyle \frac{p_i^2}{2\mu }}{\displaystyle \frac{\lambda _L}{4p_i^2}}{\displaystyle \underset{k}{}}Q_0^{}(y1)\left(1+p_k^2\right)wt_k,`$ (77)
$`{\displaystyle \frac{\lambda _L}{4p_i}}{\displaystyle \frac{l(l+1)}{2}}{\displaystyle \underset{k}{}}{\displaystyle \frac{Q_0(y1)}{p_k}}\left(1+p_k^2\right)wt_k+{\displaystyle \frac{\lambda _L\pi }{4p_i}}l(l+1)`$
$`{\displaystyle \frac{\lambda _L}{4p_i^2}}w_l^{}(1)\left(1+p_i^2\right)wt_i`$
$`K_{ij}`$ $`=`$ $`{\displaystyle \frac{\lambda _L}{4p_ip_j}}Q_l^{}(y)\sqrt{\left(1+p_i^2\right)\left(1+p_j^2\right)wt_iwt_j},(ij).`$ (78)
Despite our method to control numerical noise, numerical errors still manifest themselves in the form of spurious large negative eigenvalues for $`l8`$. Fortunately the rest of the positive eigenvalues are accurate. Some sample eigenvalues for $`0l5`$ are shown in Table \[III\] which also compares the eigenvalues generated by both the $`p`$-space and $`r`$-space codes. The $`r`$-space eigenvalues are calculated by solving NRSE using the relaxation method .

The Lande-subtraced momentum space NRSE equation with a Coulomb potential is also given in reference  and is simplified as

$`{\displaystyle \frac{p^2}{2\mu }}\varphi _{nl}(p)+{\displaystyle \frac{\lambda _C}{\pi p}}{\displaystyle _0^{\mathrm{}}}Q_l(y)\varphi _{nl}(p^{})p^{}𝑑p^{}`$ (79)
$`{\displaystyle \frac{\lambda _C}{\pi }}p\varphi _{nl}(p){\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{Q_0(y)}{p^{}}}𝑑p^{}+{\displaystyle \frac{\lambda _C\pi }{2}}p\varphi _{nl}(p)=E_{nl}\varphi _{nl}(p).`$ (80)
The kernel of a Coulomb potential can be symmetrized the same way as that of a linear potential. The matrix elements are

$`K_{ii}`$ $`=`$ $`{\displaystyle \frac{p_i^2}{2\mu }}{\displaystyle \frac{\lambda _C}{4}}p_i{\displaystyle \underset{k}{}}{\displaystyle \frac{Q_0(y1)}{p_k}}\left(1+p_k^2\right)wt_k+{\displaystyle \frac{\lambda _C\pi p_i}{2}}`$ (82)
$`{\displaystyle \frac{\lambda _C}{4}}w_l(1)\left(1+p_i^2\right)wt_i`$
$`K_{ij}`$ $`=`$ $`{\displaystyle \frac{\lambda _C}{4}}Q_l(y)\sqrt{\left(1+p_i^2\right)\left(1+p_j^2\right)wt_iwt_j},(ij).`$ (83)
The correction method which we have developed for the linear potential cannot be used in the Coulomb case. The only available technique of refining the eigenvalues of a Coulomb potential is by the way of increasing the number of mesh steps $`N`$. Some sample eigenvalues are shown in Table \[IV\]. Since both the linear and Coulomb potentials can be symmetrized using the same formalism, we can easily splice the two kernels together to calculate the eigenvalues of the Cornell (linear plus Coulomb) potential

$$V(r)=\frac{\lambda _C}{r}+\lambda _Lr.$$
(84)
It is not surprising that the correction method derived for the linear potential may also work for the Cornell potential because we expect that the error of the Cornell potential is dominated by the error of the linear potential term. But it is a surprise that the correction method works more accurately with the Cornell potential than the linear potential as evidenced by vanishingly small variance and $`\chi ^2`$.

## VIII Conclusion

The basis function method requires a priori knowledge of the eigenfunctions in order to pick out an appropriate set of basis functions. The advantage of the Nystrom method is that no such prior knowledge of the eigenfunctions is needed. The kernel constructed by the Nystrom method is also much simpler than that by the basis function method. The eigenfunctions can be generated by the same Nystrom routines that compute the eigenvalues. The Nystrom plus correction is more accurate than the basis function method in the cases studied in this paper. In other words, the new method has all of the advantages–elegance, accuracy and versatility. In addition, the kernel of the relativisitc and non-relativistic equation of motion with the Coulomb and linear potential can be symmetrized in exactly the same manner. It allows the calculation of the eigenvalues of a Cornell potential readily. Since the Nystrom method can be generalized for higher $`l`$’s, we can use it to calculate the Regge trajectories. Since the Thompson equation which we have solved is a 2-body equation, we can use it to analyze the experimental meson Regge trajectories . This will be pursued in later work.

## ACKNOWLEDGMENTS

We thank Prof. Khin Maung Maung for his helpful comments and George Nill and Daniel R. Shillinglaw for their participation.


text:
# FIRBACK II. Data Reduction and Calibration of the 170 𝜇𝑚 ISOBased on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA. Deep Cosmological Survey

## 1 Introduction

With more than 150 hours of observations, FIRBACK (Far-InfraRed BACKground) is one of the largest observational programmes made with the ISOPHOT instrument on board the Infrared Space Observatory (ISO) satellite. This cosmological deep survey at 170 $`\mu `$m covers more than 4 square degrees located in two Northern and one Southern fields (Lagache 1998, Dole 2000a). There are 106 sources detected above the sensitivity limit (180 mJy, 4$`\sigma `$). The number of sources detected above 135 mJy ($`3\sigma `$) is 196. The first result of this survey is the high number of sources observed when compared to no, or moderate, evolution models for infrared galaxies (e.g. Dole 2000a). Preliminary results in the so-called Marano1 field were published in Puget et al. (1999), and on the whole survey by Dole et al. (2000b). FIRBACK also allows for the first time the detection of the CFIRB fluctuations (Lagache & Puget 2000; Puget & Lagache, 2001).

After a presentation of the observational issues of the FIRBACK survey, we present in this paper the final stage of the data reduction and calibration. We use the Phot Interactive Analysis (PIA, Gabriel et al. 1997) for the standard reduction and calibration (Sect. 3.1) with some extra developments: flat-fielding, transient corrections and reprojection (Sect. 3.2-3.7). In Sect. 4, we review our final map calibration and compare it with absolute photometric calibration from other instruments. We then conclude in Sect. 5.

## 2 FIRBACK observations

FIRBACK covers about 4 square degrees in 3 high galactic latitude fields (chosen to have very low HI column density, typically $`N_{HI}10^{20}cm^2`$), the so-called N1 (FN1), N2 (FN2) and Marano (FSM) fields. The observations were performed at an effective wavelength of 170$`\mu \mathrm{m}`$ with the C200 array (2x2 pixels) of ISOPHOT. The pixel field of view is 1.5 arcmin. We use the AOT (Astronomical Observation Template) PHT22 in the multi-pointing staring raster mode. Each FIRBACK field (Table 1) consists of several rasters of 17$`\times `$17 pixels111Except for the FSM1 field (Puget et al. 1999) where the rasters are 19$`\times `$19 pixels:

* Eleven rasters were done in the so-called FN1 field. These rasters were re-observed with a shift in position that corresponds to a fraction of an ISOPHOT pixel to provide proper sampling where possible. Observations were done from the ISO revolution 753 to 774. This field covers the same area as the ELAIS N1 15 $`\mu `$m and 90 $`\mu `$m observations (Oliver et al. 2000).
* Nine rasters were done in the so-called FN2 field. These rasters were again re-observed with a shift in position. The field covers the same area as the ELAIS N2 15 $`\mu `$m and 90 $`\mu `$m observations. These data are FIRBACK/ELAIS observations. They were done from the ISO revolution 785 to 798.
* One raster was observed in the so-called FSM1 field (formerly called Marano 1 field) during the revolution 593 and three rasters in the so-called FSM234 field (from revolution 739 to 744). Rasters have been observed four times. Displacements between the four independent observations correspond to about 0.5 pixel, except for the FSM1 field where the displacements correspond to 2 pixels (note that the overlapped surface with the “original” Marano field is less than 15$`\%`$).

Each raster is performed in the spacecraft (Y, Z) coordinate system which is parallel to the edges of the detector array, with one pixel overlap in both Y and Z direction. The exposure time is 16s per pixel and thus 128, 128 and 256s per sky position in the FN1, FN2 and in the FSM fields respectively.

We also have a PHT25 measurement in the FSM1 field. This AOT is the absolute photometry mode for PHT-C, in which photometric calibration is achieved by chopping against the internal fine calibration source. For low fluxes (i.e. in our case), chopping is also done against the switched-off “Fine Calibration Source”, which has a temperature level of about 4K. This is the temperature of the optical support on which the instrument is mounted. Such a measurement serves to define a zero point. This mode is especially well suited for observations aiming at accurate determination of the absolute brightness of the background emission.

## 3 Data reduction

### 3.1 PIA data reduction and calibration

We use PIA, the ISOPHOT Interactive Analysis software version 7.2.2 (Gabriel et al. 1997), to correct for instrumental effects, the glitches induced by cosmic particles and to provide an initial calibration. First we apply the non-linearity correction due to 2 independent effects: the non-linearity of the Cold Readout Electronics and the downward curving ramp due to de-biasing. Deglitching is performed for each individual ramp and then the mean signal per position is derived by averaging the ramp slopes (we do not apply any transient correction at this stage). We then apply a second deglitching (for every chopper position, a “running mean” method is applied; signals that are far away by a number of sigma, typically 5, are flagged). The dark current, which represents less than 5$`\%`$ of the signal, is subtracted using the orbit dependent calibration files. We also apply the “reset interval correction”, a correction which represents less than 6$`\%`$ of the signal.

We do exactly the same data reduction for the two “Fine Calibration Source” (FCS) measurements. The FCS measurements before and after the observation vary by less than 3$`\%`$ for most of the rasters. We decide to perform the first calibration (from V/s to MJy/sr) by deriving the mean value of the two FCS measurements. We prefer not to use the interpolation between the two FCSs when they are different since it may induce Long Term Transients (LTT) that are not necessary real. On the contrary, when real, the interpolation bewteen the two FCSs could correct for the LTT, but we prefer to correct the LTT following Sect.3.3.

At this stage of the data reduction, we have the signal variation (in MJy/sr) as a function of time (or position) for each pixel in each raster. We see from Fig. 1 that the instrumental noise is very low. We have, per sky position, a very high reproductibility of the data.

### 3.2 The PHT25 measurement

The AOT PHT25 absolute measurement in the FSM1 field has the following characteristics (Fig. 2):

* 256s for the “sky dark”: observation of the cold opaque filter
* 128s for the “FCS dark”: observation of the cold FCS (turned off)
* 256s on sky
* 256s on the FCS

With this measurement we check that:

* The FCS value after 256s and the one obtained with 32s integration time (which is the integration time in the PHT22 observations) differs by less than 4$`\%`$
* The difference between the observed dark current (“sky dark”) and the default PIA one is less than 10$`\%`$
* The difference between the FCS dark current (“FCS dark”) and the default PIA one is less than 2$`\%`$

We obtain for the sky measurement a brightness222We correct the brightness for the flat-field using the average values given in Sect. 3.4 and apply the extended calibration factor correction (Sect 4.3). of 3.10$`\pm `$0.14 MJy/sr, to be compared to the value measured on the final map (Fig. 6), which is 3.08 MJy/sr. We have a very good agreement between both measurements; we do not observe any significant discrepancy between the “raster mode” observations and the “absolute mode” one.

### 3.3 Long term transients

Long Term Transients (LTT) still have an unknown origin. They are observed on a characteristic time scale of about 60-80 minutes and are the main sensitivity limitation of the ISOCAM diffuse emission analysis (Miville-Deschênes et al. 2000). LTT are observed in FIRBACK data (Fig. 3) but only for some rasters and pixels. We first tried to correct the LTT using the method developped by Miville-Deschêne et al. (2000) for ISOCAM; however, this method fails since (1) ISOCAM and ISOPHOT detectors are different and (2) the redundancy in our observations is too low. The method developed for ISOCAM clearly detects the long term variation. However, for ISOPHOT, all detectors are not contaminated by the LTT (we see Fig. 3 that the LTT is not seen by the four pixels), which is not the case for ISOCAM, where all detectors have at a first order the same behaviour. Therefore, we correct the LTT by assuming that the pixel, in each raster, that exhibits the smallest long term variation (or smallest slope) is representative of the sky. Fig. 3 shows the results before and after the correction.

### 3.4 Individual raster flat-fielding

The flat-field correction, for each individual raster, is directly estimated using the redundancy of our observations and comparing position by position the brightness of the 4 pixels. First, for each raster, we compute the responsivity of the four pixels with respect to the mean responsivity (Fig. 4) for each sky position. We see from Fig. 4 that each pixel has some high frequency instrumental noise. What we are interested in is the low-frequency difference of responsivity between the four pixels. Therefore, we derive the flat-field correction by computing a 30-points running median (and the mean of the correction is equal to 1 per sky position).

The correction is rather constant for the whole 56 rasters and equal on average to: $`1.04\pm 0.02`$, $`0.91\pm 0.02`$, $`1.09\pm 0.02`$ and $`0.94\pm 0.02`$ for pixel 1, 2, 3 and 4 respectively. This highly reproductable behaviour may be used to correct for the flat-fielding in any other C-160 observation that has no redundancy.

### 3.5 Transients induced by energetic cosmic rays

The procedure described in Sect. 3.4 also reveals some long term transients induced by energetic cosmic rays, as can been seen in Fig. 4 for the dashed line signal around the sky position number 130. These effects induced by energetic cosmic rays are quite rare: we detect only 13 events for FIRBACK (to be compared to 4 detectors $`\times `$ 56 rasters = 224 timelines). After such a cosmic ray, the signal may show an exponential-like behaviour (with a variable time constant). We correct for such events manually before computing the flat-field correction.

### 3.6 Short term transients

After a sudden increase of flux, most infrared detector signals shows an instantaneous jump followed by a slow rise to the stabilisation level; this is the short term transient. In most observations, the stabilisation level is never reached, due to the limited integration time per chopper plateau.

For the ISOPHOT C200 camera, we have information on the short term transients using our PHT25 absolute measurements. We see in Fig. 5 the signal behaviour of the PHT25 sky measurements for the 4 pixels. These measurements follow the “FCS dark” observation, whose signal is 10 times weaker (see Fig. 2). We clearly see in Fig. 5 a short term transient before the stabilisation. The instantaneous signal jump is of the order of 85$`\%`$; and the stabilisation is reached in several tens of seconds.

For FIRBACK field observations, the short term transient effect is limited due to the flatness of the fields. Flux variations are very smooth and weak, except when point sources are observed. Therefore, no correction is applied for the diffuse emission. In return, for point sources, a correction of 10$`\%`$ is applied (this 10$`\%`$ is derived using the flux measured after 16s FIRBACK integration time compared to the stabilised flux).

For the FCS FIRBACK measurements, no correction is applied since (1) the difference between the stabilised signal in the PHT25 FCS measurement and the one measured in 32s (which is the integration time on the FCS) is less than 4 $`\%`$ and (2) the second FCS observation is at the same level as the FIRBACK sky observation (no flux step) and the two FCS measurements (before and after one observation) differ by less than 3$`\%`$ for most of the rasters.

### 3.7 General “flat-fielding”

The different rasters in the same field have been done at different time (there can be several months between the rasters). Therefore, the absolute calibration performed using the FCS measurements can be slightly different from one raster to another. This difference is observed in the final FIRBACK maps as discontinuities between individual rasters.

To correct for this small absolute calibration difference, we apply a “general flat-field” correction in order to make consistent the level of each raster with its neighbourhood. This correction consists of adding or subtracting offsets, the sum of the offsets being 0. This effect is very small, less than 5$`\%`$, except for one raster in FN1 which shows a discontinuity with its neighbourhood of about 11$`\%`$.

### 3.8 Reprojection

The flux is finally projected on a 10”x10” coordinate grid using our own projection procedure. The procedure is the following:

* We compute each individual pixel coordinates.
* We create individual maps for each pixel and raster with 2D $`\frac{sinx}{x}`$ interpolation for the oversampling signal and bilinear interpolation for the coordinates.
* We compute for each field the final coordinate grid (10” sampling) using the individual raster coordinates
* At each final grid point, we add the observed signal and compute the weight and variance.

We have checked the accuracy of our method using simulated sky maps, containing both extended emission and point sources. Our reprojection does not change the background brightness and the sources flux (at the 1$`\%`$ level). One has, however, to note that this method is particularly well adapted to cosmological observations: weak sources on a rather flat background (for fields with very bright sources, this method may not preserve the point source photometry).

The final maps are shown in Fig. 7, 7 and 8 for the three fields respectively.

## 4  Calibration

The extended emission calibration has changed by a factor of about 2 between PIA version 6 and PIA version 7. A large part of this factor comes from the different footprint solid angle values used in the different PIA versions.

### 4.1  Footprint observations around Saturn

During the revolution 409, ISOPHOT made several tens of pointing around Saturn at distances between 4.2 and 45.4 arcmin in order to map the extended wings of the PHOT footprint at 170 $`\mu \mathrm{m}`$ during about 5 hours. Observations were made in the Y and Z satellite axis directions. Fig 9 shows the observed pattern on an IRAS 100 $`\mu `$m map. All directions were observed twice, back and forth. Observations were performed with the following AOTs: PHT 37-38-39 (sparse maps; PHT37: FCS; PHT38: sky; PHT39: FCS) for distances between 4.2 and 27.6 arcmin, and PHT25 (absolute photometry) at the largest distance. Integration times are around 150 sec for PHT38, 200 sec for PHT37-39 and 400 sec for PHT25.

We use PIA V7.2.2 for the data reduction and calibration. For each detector and each position, we keep only the second half of the data corresponding to the stabilised signal. These observations have no transients induced by cosmic rays and we do not apply any flat-field correction. Sixty seven files were used, the others being either unreadable (3 files) or saturated (3 files). Results are presented in Fig. 10 and 11 for all data in the Y and Z direction respectively. Each pixel is plotted and is used as an independent measurement.

### 4.2 Comparison of the footprint model with Saturn observations

A model for the ISOPHOT footprint has been developed at Heidelberg (Klaas et al., private communication) based on the ISOCAM footprint routines. This model includes the optical characteristics of the telescope, the primary and secondary mirrors, and the filters and detectors. We have used this program to compute the ISOPHOT footprint up to 20 arcmin with the 170 $`\mu `$m bandpass filter. Results of the model are shown in Fig. 12 together with the Saturn measurements for the Z axis. To make this comparison, we have assumed for Saturn a flux of 32000 Jy and removed the background using the PHT25 measurement.

We see a very good agreement between the model and the measurements. However, around 4.5 arcmin from Saturn, some data, coming from one Y direction scan (only one position), have a significantly higher flux than the model prediction. Data at similar distances from Saturn in other scans cannot be used due to saturation problems. Therefore, we interpret this discrepancy as due to detector saturation problems (one can note, however, that the contribution of these data points to the solid angle is lower than a few percent).

In conclusion, the ISOPHOT footprint measurements and model at 170 $`\mu `$m are in very good agreement. We therefore use, in the following, the model (from Klaas et al., private communication) as the definitive footprint at 170 $`\mu \mathrm{m}`$.

### 4.3 The extended brightness photometric correction factor

In PIA V6.5 the solid angle used to convert the flux in brightness was that of the pixel (1.88 10-7 sr); in PIA V7.2.2, it is the footprint model’s one but truncated at 4.1 arcmin (2.73 10-7 sr). The full solid angle of the footprint is equal to 3.06 10-7 sr. Therefore, there is a photometric correction factor to be applied to the extended emission:

$$B_\nu =B_\nu (piaV72)\times \frac{\mathrm{\Omega }_{piaV72}}{\mathrm{\Omega }_{footprint}}$$

$$B_\nu =B_\nu (piaV72)\times 0.89$$

After correction of the absolute calibration of the extended emission, point source fluxes (given in Dole et al. 2001) are computed using the footprint model (Klaas et al., private communication).

### 4.4 The rejection rate

Unique measurements have been done by the ISOPHOT team during the eclipse of the sun by the earth (Kranz et al. 1998; Klaas et al. 1998a; Lemke et al. 1998) by pointing at a sky region at 60o from the sun before, during, and after the eclipse. These measures reveal no signal variation, leading to an upper limit of 10-13 at 60o for the straylight rejection rate. This exceptional measurement clearly shows that there is no contribution to the flux from the far-side lobes. This demonstrates that ISO is able to make absolute measurements of the extended emission and gives a high degree of confidence to our absolute photometric calibration.

### 4.5 Comparison with DIRBE extended brightness

We can now compare the extended ISOPHOT FIRBACK brightness with the predicted one using DIRBE and HI measurements. Because of the size of the DIRBE beam, this comparison is only feasible due to the flatness of the FIRBACK fields; this flatness is observed with DIRBE on several degrees around each FIRBACK field.

The sky measurement is the sum of the zodiacal light, Cosmic Infrared Background (CIB) and dust interstellar emission. For each field the zodical light, at the time of the observation, has been computed using the Reach et al. (1995) DIRBE model. The CIB is extrapolated at 170 $`\mu `$m using the DIRBE measurements of Lagache et al. (2000). For the dust emission, we compute its contribution using (1) the HI column density from the Leiden-Dwingeloo survey, Hartmann & Burton 1997 (we prefer to use the HI column density rather than the DIRBE brightness since DIRBE data are very noisy in FIRBACK fields) and (2) the emissivity from Lagache et al. (2000). The final predicted emission at 170 $`\mu `$m for the three fields is shown in Table 2. It is in remarkable agreement with the measured ISOPHOT brightness333The difference between Puget et al. 2000 (paper I) FSM1 brightness and the brightness obtained now is less than 6$`\%`$.

### 4.6 Comparison with IRAS point sources measurements

The absolute flux calibration of ISOPHOT is derived using calibration standards such as planets, asteroids or stars (Klaas et al. 1998b; Schulz et al. 1999; ISO Consortium 2000a, 2000b). One can however check, using IRAS detected sources in FIRBACK fields, the consistency between IRAS 60 and 100 $`\mu `$m flux and the ISOPHOT 170 $`\mu `$m one. However, we have to keep in mind that it is very difficult to extrapolate the IRAS fluxes at 170 $`\mu `$m (one can have a factor 2 to 3 in the extrapolation using different models or black-body temperatures).

Twelve FIRBACK sources, well identified as very nearby non interacting galaxies, have IRAS 60 and 100 $`\mu `$m flux counterparts, measured with SCNAPI444SCANPI is described on the WEB site $`http://www.ipac.caltech.edu/ipac/iras/scanpi`$_$`over.html`$. We extrapolate the 60 and 100 $`\mu `$m IRAS flux using the template spectra of Dale et al. (2000). These spectra are based on the IRAS 60/100 color ratio and are thus well adapted. The comparison is shown on Fig. 13. Uncertainties on the extrapolation, that take into account only the errors on the 100 and 60 $`\mu `$m fluxes (and not the model uncertainty), are very large. Data are compatible with a slope of unity; but this comparison is only illustrative.

## 5 Conlusions

We have presented in this paper the final FIRBACK data reduction using: (1) The Phot Interactive Analysis version 7.2.2 and (2) extra developments (corrections of the flat-field, long and short term transients, transients induced by cosmic rays and adapted reprojection). Most of these extra developments have been made possible by the perfect redundancy inside each raster (one pixel overlap in both Y and Z direction). We have then checked the absolute calibration using PHT25 measurements and the footprint measured on Saturn (compared to the footprint model). We have shown that the ISOPHOT 170 $`\mu `$m extended emission calibration has to be corrected by a factor 0.89, which comes from the difference in solid angle between the PIA v7.2 and the modeled footprint.

Using this data reduction and calibration, we obtain an absolute calibration which is in remarquable agreement (better than 10$`\%`$) with brightness extrapolation that uses DIRBE data and HI column density measurements. We have in FIRBACK fields a very high signal to noise ratio (greater than 50). The main limitation comes in fact from the extragalactic source confusion itself (Dole et al. 2001).

Acknowledgements: We would like to thank A. Abergel, C. Gabriel, U. Klaas, M-A. Miville-Deschênes and the ISOPHOT team for many interactions concerning the data reduction and the footprint analysis. We thanks Danny Dale for his help in using the template spectra for IRAS extrapolation. We thank J.L. Puget for his helpful advice all along the data reduction process.