text:
# Untitled Document

CERN-TH/2003-088

Lattice QCD and the Schwarz alternating procedure

Martin Lüscher

CERN, Theory Division

CH-1211 Geneva 23, Switzerland

Abstract

A numerical simulation algorithm for lattice QCD is described, in which the short- and long-distance effects of the sea quarks are treated separately. The algorithm can be regarded, to some extent, as an implementation at the quantum level of the classical Schwarz alternating procedure for the solution of elliptic partial differential equations. No numerical tests are reported here, but theoretical arguments suggest that the algorithm should work well also at small quark masses.

1. Introduction

The simulation algorithms for (unquenched) lattice QCD that are currently in use rapidly become inefficient on large lattices and at small quark masses, where the effects of spontaneous chiral symmetry breaking set in. In the case of the HMC \[?\], the PHMC \[?,?\] and the Multiboson \[?\] algorithms, the principal technical difficulty derives from the fact that the light quark masses have to be scaled proportionally to the square of the pion mass in the chiral limit. The lattice Dirac operator is then increasingly ill-conditioned, which affects all these algorithms in a similar and rather direct way, because they all start from a global pseudo-fermion representation of the quark determinant that involves an exact (or an accurate approximate) inversion of the Dirac operator.

In the present paper a simulation algorithm is proposed that exploits the underlying local structure of the theory and that may be expected to scale in a more favourable way in the chiral regime. The general strategy is closely related to the alternating procedure that was invented by the mathematician Schwarz in the 19th century to establish the existence of the solution of the Dirichlet problem

$$\begin{array}{c}\text{ }\mathrm{\Delta }f(x)|_{x\mathrm{\Omega }}=0,f(x)|_{x\mathrm{\Omega }}=g(x),\text{(1.1)}\text{ }\hfill \end{array}$$

on arbitrary bounded domains $`\mathrm{\Omega }`$ in the plane \[?\] (for an introduction to the subject in the context of discretized partial differential equations see ref. \[?\], for example). Very briefly this method obtains the solution iteratively by dividing $`\mathrm{\Omega }`$ into a set of overlapping subdomains and by solving the Dirichlet problem on these, in each step of the iteration, with boundary values determined from the current approximation to the solution.

The design of simulation algorithms for lattice QCD that operate on overlapping blocks of lattice points is non-trivial, however, because the global correctness of the simulation must be guaranteed. In principle the problem can be solved using stochastic acceptance–rejection steps similar to those previously considered by Hasenbusch \[?\] (see also refs. \[?–?\]). The key question is then whether these correction steps can be implemented so that high acceptance rates are achieved, and a significant part of the present paper is therefore devoted to this issue.

2. General form of the algorithm

In this section the proposed algorithm is described in outline. Most technical details are deferred to the later sections and important improvements (such as preconditioning) are omitted in order to keep the presentation as simple as possible. Some algorithm theory, as summarized in appendix A, is nevertheless required to be able to understand the procedure.

2.1 Preliminaries

Although the algorithm is more generally applicable, only the standard Wilson formulation \[?\] of lattice QCD will be considered here (optionally including $`\mathrm{O}(a)`$ improvement \[?,?\]) with a doublet of mass-degenerate quarks. The lattice spacing is set to unity for convenience and the SU(3) link variables are denoted by $`U(x,\mu )`$ as usual.

After integration over the quark fields, the probability distribution of the gauge field reads

$$\begin{array}{c}\text{ }P[U]=\frac{1}{𝒵}\mathrm{e}^{S_\mathrm{G}}detQ^2,Q\gamma _5(D_\mathrm{w}+m_0),\text{(2.1)}\text{ }\hfill \end{array}$$

where $`𝒵`$ denotes the partition function, $`S_\mathrm{G}`$ the plaquette action, $`D_\mathrm{w}`$ the Wilson–Dirac operator and $`m_0`$ the bare quark mass.

2.2 Alternating procedure

Following the classical Schwarz procedure, the proposed algorithm visits rectangular blocks of lattice points according to some scheme and updates the link variables residing there. The blocks can have arbitrary sizes in principle, and their position may be chosen randomly, for example, so that all link variables are treated equally. However, as will become clear later, the algorithm is designed to perform particularly well if the blocks are small in physical units, i.e. if their edges are less than about $`1`$ fm long.

On each block $`\mathrm{\Lambda }`$ that is visited in the course of this process, changes in the gauge field on $`\mathrm{\Lambda }`$ are proposed that satisfy detailed balance with respect to the distribution

$$\begin{array}{c}\text{ }P_\mathrm{\Lambda }[U]=\frac{1}{𝒵_\mathrm{\Lambda }}\mathrm{e}^{S_\mathrm{G}}det(Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{})^2.\text{(2.2)}\text{ }\hfill \end{array}$$

A stochastic acceptance–rejection step then needs to be applied to correct for the difference between this distribution and the exact distribution (2.1). The operators $`Q_\mathrm{\Lambda }`$ and $`Q_\mathrm{\Lambda }^{}`$ that appear here coincide with $`Q`$, except that they act on Dirac fields defined on the block $`\mathrm{\Lambda }`$ and its complement $`\mathrm{\Lambda }^{}`$, respectively, with Dirichlet boundary conditions. At the level of the fermion action, a complete decoupling of the quark fields inside and outside the block is thus achieved. In particular, the proposals for the link variables residing in $`\mathrm{\Lambda }`$ can be generated locally using a block version of the HMC algorithm, for example.

The stochastic acceptance–rejection step, on the other hand, involves an inversion of the Dirac operator $`Q`$ on the full lattice. There is a fair amount of choice in the detailed implementation of this step, which can be exploited to reduce the influence on the acceptance probability of the link variables far away from the block. Moreover, the suggested procedure (which is explained in sect. 5) restricts the pseudo-fermion field that needs to be introduced at this stage to the boundary of the block. At least to some extent, the local character of the algorithm is thus preserved.

2.3 Hierarchical structure

In the form described above the algorithm requires a computational effort per block update that increases roughly proportionally to the lattice volume. By introducing a hierarchy of blocks such as the one shown in fig. 1, one may, however, be able to do better than this. Proposed changes of the link variables on the smallest blocks are then obtained as before, and these are taken as proposed configurations on the next larger blocks, and so on, accepting or rejecting them so that detailed balance is satisfied with respect to the associated distributions $`P_\mathrm{\Lambda }[U]`$. Finally a simultaneous global acceptance–rejection step is applied to the surviving configurations.

Fig. 1. Proposed updates of the link variables on small blocks (dark grey squares) can be filtered through acceptance–rejection steps on increasingly larger blocks (light grey squares) before the global acceptance–rejection step is applied.

This procedure is logically correct since non-overlapping blocks are decoupled from each other and can be updated in parallel. On the other hand, it will only work out in practice if the configurations that have been accepted on the largest blocks have a high probability to be accepted on the whole lattice. For the chosen implementation of the acceptance–rejection steps, a theoretical argument will be given in sect. 5 that explains why this should be so (under certain conditions).

2.4 Low-mode reweighting

On physically small blocks $`\mathrm{\Lambda }`$, the Dirac operator $`Q_\mathrm{\Lambda }`$ that acts on the quark fields in $`\mathrm{\Lambda }`$ is not expected to become ill-conditioned at small quark masses. The situation here is actually similar to the one encountered in studies of the Schrödinger functional, where the boundary conditions provide an infrared cutoff on both the gluon and the quark modes \[?,?\]. As a consequence the block simulation algorithm should work well even if the bare quark mass $`m_0`$ is set to the critical mass $`m_c`$ \[?\].

However, since the Wilson formulation of lattice QCD violates chiral symmetry, the Dirac operator on large lattices is not protected against the accidental presence of eigenvectors with eigenvalues much smaller than $`m_0m_c`$. When a gauge field configuration is proposed where this happens, it will only rarely pass the global acceptance–rejection step. Such configurations are therefore sampled with low statistics and this can easily lead to uncontrolled statistical fluctuations if a quantity is considered that is sensitive to the low modes of the Dirac operator (the propagator of the pseudo-scalar density, for example).

It may be possible to solve this problem by including the product

$$\begin{array}{c}\text{ }W[U]=\underset{k=1}{\overset{n}{}}\lambda _k\text{(2.3)}\text{ }\hfill \end{array}$$

of the first few eigenvalues $`\lambda _1,\mathrm{},\lambda _n`$ of $`Q^2`$ in the observables and the inverse factor $`W[U]^1`$ in the global acceptance–rejection step. The computational overhead for this modification may not be negligible, but it should be noted in this connection that the eigenvalues do not need to be computed very accurately (as long as a definite procedure is used that obtains the eigenvalues independently of the previous gauge field configurations). Moreover the associated approximate eigenvectors can help to accelerate the inversion of the Dirac operator that is required in the acceptance–rejection step \[?\].

3. Domain decomposition

Let $`\mathrm{\Lambda }`$ be an arbitrary rectangular block of lattice points. $`\mathrm{\Lambda }`$ and its complement $`\mathrm{\Lambda }^{}`$ (the set of points not in $`\mathrm{\Lambda }`$) define a particular case of a domain decomposition of the lattice. In the following paragraphs the aim is to introduce some basic notation related to this decomposition, but the terminology is quite general and extends to the case where $`\mathrm{\Lambda }`$ is replaced by the union of a set of non-overlapping blocks.

3.1 Boundary points

On the lattice it is important to distinguish between interior and exterior boundary points (see fig. 2). The boundary values for the classical Dirichlet problem on $`\mathrm{\Lambda }`$, for example, should be specified on the set $`\mathrm{\Lambda }`$ of all exterior boundary points if the standard nearest-neighbour lattice laplacian is used, while the set of interior boundary points plays an analogous rôle from the point of view of the complementary domain $`\mathrm{\Lambda }^{}`$ and is therefore denoted by $`\mathrm{\Lambda }^{}`$.

Fig. 2. Two-dimensional slice of a $`6^4`$ block (black and grey points inside the dashed square). The grey and the open points represent the interior and the exterior boundary points of the block respectively.

A special feature of this geometry is that for each exterior boundary point there is a unique closest point on the interior boundary, referred to as its partner point. The converse is evidently not true, i.e. different exterior boundary points may have the same partner point.

3.2 Decomposition of the Dirac operator

In position space the hermitian lattice Dirac operator is a sparse matrix that assumes the block-diagonal form

$$\begin{array}{c}\text{ }Q=\left(\begin{array}{cc}Q_\mathrm{\Lambda }& Q_\mathrm{\Lambda }\\ \\ Q_\mathrm{\Lambda }^{}& Q_\mathrm{\Lambda }^{}\end{array}\right)\text{(3.1)}\text{ }\hfill \end{array}$$

if the lattice points are ordered so that those in $`\mathrm{\Lambda }`$ come first. The operator $`Q_\mathrm{\Lambda }`$, for example, acts on Dirac fields on $`\mathrm{\Lambda }`$ in the same way as $`Q`$, except that all terms involving the exterior boundary points are set to zero (which is equivalent to imposing Dirichlet boundary conditions on $`\mathrm{\Lambda }`$).

It is often convenient to consider $`Q_\mathrm{\Lambda }`$, $`Q_\mathrm{\Lambda }^{}`$, $`Q_\mathrm{\Lambda }`$ and $`Q_\mathrm{\Lambda }^{}`$ to be operators in the space of Dirac fields that are defined on the whole lattice rather than on $`\mathrm{\Lambda }`$ or $`\mathrm{\Lambda }^{}`$ only. The embedding is done in the obvious way by padding with zeros. Equation (3.1) may then be written in the form

$$\begin{array}{c}\text{ }Q=Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{}+Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{},\text{(3.2)}\text{ }\hfill \end{array}$$

and an equivalent expression for the determinant in eq. (2.2) is $`detQ_\mathrm{\Lambda }^2detQ_\mathrm{\Lambda }^{}^2`$. This notation is perhaps slightly abusive but should not lead to any confusion since it is usually clear from the context which interpretation applies.

The operator $`Q_\mathrm{\Lambda }`$ and its hermitian conjugate $`Q_\mathrm{\Lambda }^{}`$ connect the domains $`\mathrm{\Lambda }`$ and $`\mathrm{\Lambda }^{}`$ to each other. Explicitly $`Q_\mathrm{\Lambda }`$ is given by

$$\begin{array}{c}\text{ }Q_\mathrm{\Lambda }\psi (x)=\gamma _5\theta _\mathrm{\Lambda }(x)\underset{\mu =0}{\overset{3}{}}\{\frac{1}{2}(1\gamma _\mu )\theta _\mathrm{\Lambda }^{}(x+\widehat{\mu })U(x,\mu )\psi (x+\widehat{\mu })\text{ }\hfill \\ \\ \text{ }+\frac{1}{2}(1+\gamma _\mu )\theta _\mathrm{\Lambda }^{}(x\widehat{\mu })U(x\widehat{\mu },\mu )^1\psi (x\widehat{\mu })\},\text{(3.3)}\text{ }\hfill \end{array}$$

where $`\theta _\mathrm{\Lambda }`$ and $`\theta _\mathrm{\Lambda }^{}`$ denote the characteristic functions of $`\mathrm{\Lambda }`$ and $`\mathrm{\Lambda }^{}`$ respectively$`\mathrm{}`$ The normalization conventions and all unexplained notations are as in ref. \[?\].. This operator thus transports the Dirac spinors on the exterior boundary points to the partner points on the interior boundary, while $`Q_\mathrm{\Lambda }^{}`$ does the same in the opposite direction.

3.3 Dirichlet problem & space of boundary values

To be able to properly pose the Dirichlet boundary value problem for the lattice Dirac operator on $`\mathrm{\Lambda }`$, the space of boundary values first needs to be specified. The situation is practically the same as in the case of the Schrödinger functional which was studied in ref. \[?\]. It is useful in this connection to introduce the projector-valued function

$$\begin{array}{c}\text{ }\theta _\mathrm{\Lambda }(x)=\underset{\mu =0}{\overset{3}{}}\theta _\mathrm{\Lambda }^{}(x)\left\{\frac{1}{2}(1\gamma _\mu )\theta _\mathrm{\Lambda }(x\widehat{\mu })+\frac{1}{2}(1+\gamma _\mu )\theta _\mathrm{\Lambda }(x+\widehat{\mu })\right\}\text{(3.4)}\text{ }\hfill \end{array}$$

that will, in many respects, play the rôle of the characteristic function of the exterior boundary. The only non-vanishing terms on the right-hand side of eq. (3.4) are in fact those where $`x`$ is in $`\mathrm{\Lambda }`$ and $`x\widehat{\mu }`$ or $`x+\widehat{\mu }`$ its partner point, and the projectors $`\frac{1}{2}(1\pm \gamma _\mu )`$ are precisely such that the identity

$$\begin{array}{c}\text{ }Q_\mathrm{\Lambda }=Q_\mathrm{\Lambda }\theta _\mathrm{\Lambda }\text{(3.5)}\text{ }\hfill \end{array}$$

holds.

The Dirichlet boundary value problem for the lattice Dirac operator is now to find a Dirac field $`\psi (x)`$ on $`\mathrm{\Lambda }\mathrm{\Lambda }`$ that satisfies

$$\begin{array}{c}\text{ }Q\psi (x)=0\text{for all }x\mathrm{\Lambda },\text{(3.6)}\text{ }\hfill \\ \\ \text{ }\theta _\mathrm{\Lambda }(x)\psi (x)=\psi (x)=\eta (x)\text{for all }x\mathrm{\Lambda },\text{(3.7)}\text{ }\hfill \end{array}$$

where $`\eta =\theta _\mathrm{\Lambda }\eta `$ is any prescribed field on $`\mathrm{\Lambda }`$. Recalling the decomposition (3.1) of the Dirac operator, it is clear that $`Q`$ can be replaced by $`Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }`$ in eq. (3.6). At all points in $`\mathrm{\Lambda }`$ the solution is then given by

$$\begin{array}{c}\text{ }\psi =Q_\mathrm{\Lambda }^1Q_\mathrm{\Lambda }\eta ,\text{(3.8)}\text{ }\hfill \end{array}$$

where use has been made of the fact that $`Q_\mathrm{\Lambda }`$ moves the field $`\eta `$ from the exterior to the interior boundary of the block. In particular, the solution is uniquely determined by the specified boundary values (if $`Q_\mathrm{\Lambda }`$ is invertible).

4. Block simulation algorithm

For the generation of proposed gauge field configurations on a given block $`\mathrm{\Lambda }`$, the HMC and PHMC algorithms are both equally suitable. Some relevant details are given in this section for the case of the HMC algorithm. The correctness of the local form of this algorithm is easily shown by noting that it matches the general pattern outlined in subsect. A.4.

4.1 Active and spectator link variables

As already mentioned in sect. 2, the proposed updates of the gauge field on $`\mathrm{\Lambda }`$ must satisfy detailed balance with respect to the distribution (2.2), which is now written in the form

$$\begin{array}{c}\text{ }P_\mathrm{\Lambda }[U]=\frac{1}{𝒵_\mathrm{\Lambda }}\mathrm{e}^{S_\mathrm{G}}detQ_\mathrm{\Lambda }^2detQ_\mathrm{\Lambda }^{}^2.\text{(4.1)}\text{ }\hfill \end{array}$$

To achieve a complete decoupling from the surrounding lattice, only the link variables residing on a subset of links in $`\mathrm{\Lambda }`$, the active link variables, should be changed. In particular, if the links shown in fig. 3 are selected, the factor $`detQ_\mathrm{\Lambda }^{}`$ is guaranteed to be independent of the active link variables, independently of whether the $`\mathrm{O}(a)`$-improved or the unimproved theory is considered.

Fig. 3. The block update algorithm only changes the link variables on a subset of links in the chosen block (solid lines). In particular, the link variables on the interior boundary and those connecting the interior to the exterior boundary are held fixed.

At this point the HMC algorithm can be set up as usual by introducing a pseudo-fermion field $`\varphi (x)`$ on $`\mathrm{\Lambda }`$ and the canonical momenta $`\mathrm{\Pi }(x,\mu )`$ of the active link variables $`U(x,\mu )`$. The total action of the system is then taken to be

$$\begin{array}{c}\text{ }S=S_\mathrm{G}+\frac{1}{2}(\mathrm{\Pi },\mathrm{\Pi })+(\varphi ,Q_\mathrm{\Lambda }^2\varphi ),\text{(4.2)}\text{ }\hfill \end{array}$$

where the brackets $`(,)`$ denote the natural scalar products in the relevant spaces of fields. Note that the factor $`detQ_\mathrm{\Lambda }^{}`$ can be ignored, since it does not depend on the active link variables. All inactive link variables actually play a spectator rôle in the following, which is rather similar to what happens in the case of the widely used simulation algorithms for pure gauge theories, where the link variables are updated one after the other.

4.2 Variable-speed molecular dynamics

Following the standard procedure, a new configuration of the active link variables is now generated by first choosing the pseudo-fermion field and the canonical momenta randomly with conditional probability proportional to $`\mathrm{e}^S`$. The momenta and the active link variables are then evolved from their initial values $`\mathrm{\Pi },U`$ to some final values $`\mathrm{\Pi }^{},U^{}`$ by solving the molecular dynamics equations

$$\begin{array}{c}\text{ }\frac{\mathrm{d}}{\mathrm{d}\tau }\mathrm{\Pi }(x,\mu )=\gamma (x,\mu )F(x,\mu ),\text{(4.3)}\text{ }\hfill \\ \\ \text{ }\frac{\mathrm{d}}{\mathrm{d}\tau }U(x,\mu )=\gamma (x,\mu )\mathrm{\Pi }(x,\mu )U(x,\mu ),\text{(4.4)}\text{ }\hfill \end{array}$$

from $`\tau =0`$ to $`\tau =1`$. As usual the force $`F(x,\mu )`$ takes values in the Lie algebra of $`\mathrm{SU}(3)`$ and is determined by the requirement that

$$\begin{array}{c}\text{ }(\omega ,F)=\underset{ϵ0}{lim}\left\{\frac{1}{ϵ}\left(S|_{UU_ϵ}S\right)\right\},U_ϵ(x,\mu )\mathrm{e}^{ϵ\omega (x,\mu )}U(x,\mu ),\text{(4.5)}\text{ }\hfill \end{array}$$

for all variations $`\omega (x,\mu )`$ of the active link variables. It should be emphasized that the computation of the force is an operation that is local to the block $`\mathrm{\Lambda }`$. In particular, only the Dirac operator $`Q_\mathrm{\Lambda }`$ needs to be inverted.

The action and the integration measure are both preserved in the course of the molecular dynamics evolution, for any choice of the (real) speed factor $`\gamma (x,\mu )`$. This factor is usually set to $`1`$, but for reasons explained later, it may be advisable, in the present context, to take a smaller value on the links close to the boundary of $`\mathrm{\Lambda }`$. The associated link variables are then slowed down and tend to deviate less from their initial values at the end of the evolution (stochastic equations where different field components are evolved at different speeds have previously been considered by Davies et al. \[?\]).

It is, incidentally, straightforward to incorporate this modification in the numerical integration scheme that is used to solve the molecular dynamics equations. As in the case of the standard HMC algorithm, an acceptance–rejection step is finally required to correct for the integration error and thus to obtain a transition probability that satisfies detailed balance with respect to the distribution $`P_\mathrm{\Lambda }[U]`$.

5. Global acceptance–rejection step

The gauge field configurations that are generated by applying the local algorithm described above can be taken as proposed configurations for the theory on a larger block or on the whole lattice. Whether a proposed field is accepted or not is then decided by applying a stochastic acceptance–rejection step. A particular implementation of this step is now going to be discussed, first for the case of the transition from a given block $`\mathrm{\Lambda }`$ to the full lattice.

5.1 Acceptance probability

Following the general strategy summarized in subsect. A.5, the probability distribution $`P_\mathrm{\Lambda }[U]`$ is considered to be an approximation to the exact distribution $`P[U]`$. An auxiliary pseudo-fermion field $`\chi `$ is then added to the system with action

$$\begin{array}{c}\text{ }S_\chi =(\chi ,(Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{})Q^1\varrho _\mathrm{\Lambda }^2Q^1(Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{})\chi ),\text{(5.1)}\text{ }\hfill \end{array}$$

where a local shape function

$$\begin{array}{c}\text{ }\varrho _\mathrm{\Lambda }(x)=\theta _\mathrm{\Lambda }(x)+ϵ\left(1\theta _\mathrm{\Lambda }(x)\right),\text{(5.2)}\text{ }\hfill \end{array}$$

with parameter $`ϵ>0`$, has been included for reasons that will become clear in a moment. The simultaneous probability distribution of the enlarged theory,

$$\begin{array}{c}\text{ }\widehat{P}[\chi ,U]=\frac{1}{\widehat{𝒵}}\mathrm{e}^{S_\chi }P_\mathrm{\Lambda }[U],\text{(5.3)}\text{ }\hfill \end{array}$$

reproduces the correct distribution when the auxiliary field is integrated out. Note that $`\varrho _\mathrm{\Lambda }`$ is invertible and independent of the fields. The associated determinant is therefore a constant that cancels out once the normalization factors are taken into account.

Before the acceptance–rejection step can be applied, a pseudo-fermion field $`\chi `$ must be generated randomly with conditional probability $`\widehat{P}[\chi |U]`$. This is easily achieved by setting

$$\begin{array}{c}\text{ }\chi =\left(Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{}\right)^1Q\varrho _\mathrm{\Lambda }^1\eta ,\text{(5.4)}\text{ }\hfill \end{array}$$

where $`\eta `$ is chosen with probability proportional to $`\mathrm{e}^{(\eta ,\eta )}`$. The proposed field $`U^{}`$ generated by the block simulation algorithm is then accepted with probability

$$\begin{array}{c}\text{ }P_{\mathrm{acc}}=\mathrm{min}\{1,\mathrm{e}^{\mathrm{\Delta }S_\chi }\},\mathrm{\Delta }S_\chi S_\chi |_{UU^{}}S_\chi ,\text{(5.5)}\text{ }\hfill \end{array}$$

which ensures that the gauge field is correctly updated with a transition probability that satisfies detailed balance with respect to the exact distribution $`P[U]`$.

5.2 Elimination of the volume terms

In the following lines it is shown that $`\mathrm{\Delta }S_\chi `$ only depends on the values of $`\eta `$ on the exterior boundary $`\mathrm{\Lambda }`$, up to terms that are proportional to the parameter $`ϵ`$ of the shape function $`\varrho _\mathrm{\Lambda }`$ and that can, therefore, be eliminated by taking the limit $`ϵ0`$.

The proof starts from the identities

$$\begin{array}{c}\text{ }\mathrm{\Delta }S_\chi =(\zeta ,\zeta )(\eta ,\eta ),\text{(5.6)}\text{ }\hfill \\ \\ \text{ }\zeta \varrho _\mathrm{\Lambda }Q^1\left(Q_\mathrm{\Lambda }^{}+Q_\mathrm{\Lambda }^{}^{}\right)\left(Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{}\right)^1Q\varrho _\mathrm{\Lambda }^1\eta \text{ }\hfill \\ \\ \text{ }=\eta +\varrho _\mathrm{\Lambda }Q^1\left(Q_\mathrm{\Lambda }^{}Q_\mathrm{\Lambda }\right)\left(Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{}\right)^1\left(Q_\mathrm{\Lambda }+Q_\mathrm{\Lambda }^{}\right)\varrho _\mathrm{\Lambda }^1\eta ,\text{(5.7)}\text{ }\hfill \end{array}$$

where the last line has been obtained using the block decomposition (3.2) of the Dirac operator and the fact that $`Q_\mathrm{\Lambda }^{}`$, $`Q_\mathrm{\Lambda }`$ and $`Q_\mathrm{\Lambda }^{}`$ are independent of the active link variables. The crucial observation is then that all terms proportional to $`Q_\mathrm{\Lambda }^{}\varrho _\mathrm{\Lambda }^1\eta `$ vanish since $`Q_\mathrm{\Lambda }^{}Q_\mathrm{\Lambda }`$ acts in the subspace of fermion fields supported on $`\mathrm{\Lambda }`$. Recalling eqs. (3.5) and (5.2), this leads to

$$\begin{array}{c}\text{ }\zeta =\eta +\varrho _\mathrm{\Lambda }Q^1\left(Q_\mathrm{\Lambda }^{}Q_\mathrm{\Lambda }\right)Q_\mathrm{\Lambda }^1Q_\mathrm{\Lambda }\eta \text{(5.8)}\text{ }\hfill \end{array}$$

and thus to an expression

$$\begin{array}{c}\text{ }\mathrm{\Delta }S_\chi =2\mathrm{Re}(\theta _\mathrm{\Lambda }\eta ,\xi )+(\xi ,\xi )+\mathrm{O}(ϵ),\text{(5.9)}\text{ }\hfill \\ \\ \text{ }\xi \theta _\mathrm{\Lambda }Q^1\left(Q_\mathrm{\Lambda }^{}Q_\mathrm{\Lambda }\right)Q_\mathrm{\Lambda }^1Q_\mathrm{\Lambda }\eta ,\text{(5.10)}\text{ }\hfill \end{array}$$

of the type announced above.

At this point the parameter $`ϵ`$ can safely be taken to zero because all terms that are inversely proportional to $`ϵ`$ have disappeared. Moreover, since only the component $`\theta _\mathrm{\Lambda }\eta `$ of the field $`\eta `$ enters the final expressions, the constraint

$$\begin{array}{c}\text{ }\eta =\theta _\mathrm{\Lambda }\eta \text{(5.11)}\text{ }\hfill \end{array}$$

may now be imposed without loss. The field is thus restricted to the linear space of boundary values for the Dirichlet problem on $`\mathrm{\Lambda }`$ (cf. subsect. 3.3).

5.3 Acceptance rate

Evidently the whole approach will fail in practice unless the acceptance probability (5.5) is reasonably large on average. This will obviously be the case if $`\mathrm{\Delta }S_\chi `$ is only rarely greater than $`1`$, and the question thus arises of whether there is any theoretical reason to expect this to be so.

Some insight into the problem can be gained by assuming, for a moment, that $`U^{}`$ differs from $`U`$ only on a single link $`l`$ in $`\mathrm{\Lambda }`$. In eq. (5.9) the first term then involves the propagation of the random field $`\eta `$ from the exterior boundary $`\mathrm{\Lambda }`$ to the endpoints of $`l`$ and a second propagation from there back to the boundary. The other term has a similar structure but involves altogether four quark propagators.

Quark propagators typically decay like $`d^3`$ in the short-distance regime and more rapidly at larger distances $`d`$. The contributions to $`\mathrm{\Delta }S_\chi `$ with two quark propagators are thus suppressed by a factor $`d^3d^3`$ at least, where $`d`$ and $`d^{}`$ are the distances of the link $`l`$ from the initial and final points on $`\mathrm{\Lambda }`$. The sum over all these points then still needs to be performed, but since $`\eta `$ is a random field there are strong cancellations in this sum and on average the result is enhanced by a factor proportional to the number of boundary points only (rather than its square). An important suppression factor thus remains, particularly when the link $`l`$ is a fair distance away from the boundary.

Reasonable acceptance rates can thus be expected if the proposed changes in the link variables on the links closer to the boundary of the block are not too large. By adjusting the speed factor $`\gamma (x,\mu )`$ in the molecular dynamics equations (4.3),(4.4), this can be easily achieved. Note that there must also be important cancellations in the sum over the active link variables that is implicit in eq. (5.10), because they change in random directions in group space.

5.4 Numerical exercises

Some numerical confirmation of this qualitative argumentation is clearly desirable at this point. The decay of the quark propagator away from the boundary of the block, for example, can be checked by calculating the solution $`\psi =Q_\mathrm{\Lambda }^1Q_\mathrm{\Lambda }\eta `$ of the Dirichlet problem on $`\mathrm{\Lambda }`$. A representative result of such studies for the average of $`\left|\psi (x)\right|^2`$ over all boundary values and gauge fields is reported in fig. 4.

Another check on the correctness of the theoretical argumentation is obtained by computing the average of the acceptance probability $`P_{\mathrm{acc}}`$ in the quenched approximation, where the proposed changes in the active link variables are generated using a link update algorithm. The results in this case are rather encouraging too and they demonstrate, in particular, the importance of the elimination of the volume terms.

Fig. 4. Average magnitude squared of the solution $`\psi (x)`$ of the Dirichlet problem (3.6),(3.7) on a $`12^4`$ block as a function of the distance $`d`$ from the exterior boundary of the block. The data shown are from a simulation of quenched QCD at a lattice spacing of about $`0.1`$ fm. Diamonds, squares and circles correspond to values of the quark mass where the pion mass is known to be approximately equal to $`320`$, $`457`$ and $`579`$ MeV respectively \[?\].

In all these studies the dependence of the calculated quantities on the quark mass appears to be fairly weak, as long as the edges of the block are smaller than $`1`$ fm or so. This is perhaps not totally surprising in view of what has been said above, but the observation suggests that the algorithm will remain effective also at small quark masses.

5.5 Hierarchical filter

As was briefly discussed in subsect. 2.3, a hierarchical structure should be adopted on large lattices where the proposed gauge field configurations are filtered through a sequence of blocks of increasing size. The filter involves an acceptance–rejection step for each transition from a block $`\mathrm{\Lambda }`$ to the next larger block $`\mathrm{\Omega }`$. Equations (5.9)–(5.11) remain valid in this case if $`Q^{}`$ is replaced by $`Q_\mathrm{\Omega }^{}`$. An important point to note here is that the distance of the links where the active link variables reside to the boundary of the blocks is increasing. The acceptance probability for the step from one block to the next is therefore expected to rise quickly.

The proposed configurations that remain after the application of the hierarchical filter must finally be submitted to a global acceptance–rejection step. If several configurations are submitted simultaneously, as suggested in subsect. 2.3, the acceptance probability is again given by eqs. (5.5) and (5.9)–(5.11), where $`\mathrm{\Lambda }`$ now stands for the union of the big blocks that contain the proposed configurations. The boundary $`\mathrm{\Lambda }`$ is then the union of the boundaries of these blocks and $`\theta _\mathrm{\Lambda }`$ the sum of the corresponding characteristic functions. In particular, there are terms in $`\mathrm{\Delta }S_\chi `$ that couple the random field $`\eta `$ on the boundaries of different blocks.

6. Miscellaneous remarks

(1) Preconditioning. As is generally the case in lattice QCD, the preconditioning of the Dirac operator can be expected to have a significant impact on the performance of the proposed algorithm. In particular, it is straightforward to implement even–odd preconditioning both in the block simulation algorithm and in the global acceptance–rejection step.

(2) Block simulation algorithm. There is a range of algorithms that can be used to generate the proposed gauge field configurations on a given block. In particular, most improvements of the HMC and the PHMC algorithms that have been introduced over the years, such as the two-boson method \[?,?\], should be beneficial here too.

A more radical possibility is to generate the configurations through a standard link update algorithm, using an adapted gauge field action, and to apply a stochastic acceptance–rejection step to correct for the quark determinant \[?–?\]. Following the lines of subsect. 5.2, the random pseudo-fermion field that needs to be introduced in the correction step can be reduced to the support of $`Q_\mathrm{\Lambda }^{}Q_\mathrm{\Lambda }`$ in this case. This will no doubt increase the average acceptance probability, but whether the approach can compete with the block HMC algorithm, for example, remains to be seen.

(3) Parallel processing. The proposed algorithm is well suited for parallel computers, because distant blocks can be updated independently of each other. Different distributions of the workload are conceivable, and it is not required that the smallest blocks be processed on single nodes, although in this case the communication overhead is minimized. It is, incidentally, advisable to translate the gauge field after each cycle rather than shifting the block positions. The program then always operates on the same blocks and their position can be chosen so as to fit the processor grid in the best possible way.

7. Conclusions

The problem to find efficient simulation algorithms for unquenched lattice QCD in the chiral regime has been around for many years. Whether the approach described in this paper provides a viable solution is not certain and can only be decided after extensive numerical tests have been performed.

It is quite clear, however, that significant progress in this area can only be made if the short- and long-distance effects of the sea quarks are treated differently. This insight is not new and has motivated the truncated determinant approximation of Duncan et al. \[?,?\], for example, and a modification of the HMC algorithm with multiple molecular dynamics time scales \[?\]. To some extent, the two-boson method studied in refs. \[?,?\] can also be put under this general heading.

The algorithm proposed here separates the short-distance effects from the long-distance ones by updating the gauge field on physically small blocks of lattice points that are visited sequentially. On the blocks the theory can be simulated even at vanishing quark masses, because the chosen boundary conditions imply an infrared cutoff on the spectrum of the Dirac operator. The long-distance effects of the sea quarks are then incorporated by a sequence of acceptance–rejection steps that lead from the small blocks to larger blocks and eventually to the full lattice.

I wish to thank Martin Hasenbusch and Ulli Wolff for correspondence on fermion simulation algorithms and Peter Weisz for a critical reading of the manuscript. The computer time for the numerical studies reported in sect. 5 has been kindly provided by DESY–Hamburg and by the Institute for Theoretical Physics at the University of Berne.

Appendix A. Transition probabilities

To prove the correctness of the algorithm proposed in this paper it suffices to write down the transition probabilities that are associated to the various algorithmic steps and to recall some basic facts from the general theory of numerical simulations. For completeness, the relevant theoretical results are summarized in this appendix.

A.1 Convergence theorem

To avoid any unnecessary complications, an abstract discrete system will be considered with a finite number of states $`s`$. The thermal equilibrium distribution of the states is denoted by $`P(s)`$ and it is assumed that $`P(s)>0`$ for all $`s`$.

Numerical simulation algorithms for this system (as far as they follow the standard procedures) are defined by a transition probability $`T(ss^{})`$ that satisfies

| $`1.`$ | $`T(ss^{})`$ is non-negative and such that $`_s^{}T(ss^{})=1`$ for all $`s`$. |
| --- | --- |
| $`2.`$ | The equilibrium distribution is preserved, i.e. for all states $`s^{}`$ the equation $`_sP(s)T(ss^{})=P(s^{})`$ holds. |
| $`3.`$ | The recurrence probability $`T(ss)`$ is non-zero for all $`s`$ and any state can be reached from any other state in a finite number of transitions. |

In practice the simulation starts from an initial state $`s_0`$ and generates a sequence of states $`s_1,s_2,\mathrm{}`$ recursively with transition probability $`T(s_ks_{k+1})`$. The fundamental theorem of simulation theory then asserts that the states in the sequence will be asymptotically distributed according to the equilibrium distribution.

A.2 Composition rule & detailed balance

If $`T_1`$ and $`T_2`$ are any two transition probabilities that satisfy conditions 1–3, it is trivial to show that the same is true for the composed transition probability

$$\begin{array}{c}\text{ }T(ss^{})=\underset{r}{}T_1(sr)T_2(rs^{}).\text{(}\text{A}\text{.1)}\text{ }\hfill \end{array}$$

The associated algorithm generates an intermediate state $`r`$ with probability $`T_1(sr)`$ and then the new state $`s^{}`$ with probability $`T_2(rs^{})`$. Evidently any number of algorithms can be combined in this way.

Transition amplitudes that satisfy detailed balance,

$$\begin{array}{c}\text{ }P(s)T(ss^{})=P(s^{})T(s^{}s),\text{(}\text{A}\text{.2)}\text{ }\hfill \end{array}$$

are an important special case. If detailed balance holds, property 2 is an immediate consequence of property 1. However, the composed transition probability (A.1) in general does not satisfy detailed balance if $`T_1`$ and $`T_2`$ do, because the order of the factors in eq. (A.1) matters. It is possible to correct for this deficit by choosing the order randomly or by forming reversion-symmetric products.

A.3 Acceptance--rejection method

Valid algorithms may often be obtained from transition amplitudes $`T_0(ss^{})`$ that satisfy detailed balance with respect to some approximate distribution $`P_0(s)`$ \[?\]. The idea is to first propose a new state $`s^{}`$ according to this probability and to accept it with probability

$$\begin{array}{c}\text{ }P_{\mathrm{acc}}(s,s^{})=\mathrm{min}\{1,\frac{P_0(s)P(s^{})}{P(s)P_0(s^{})}\}.\text{(}\text{A}\text{.3)}\text{ }\hfill \end{array}$$

Otherwise (i.e. when $`s^{}`$ is rejected) the old state $`s`$ is taken to be the new one.

The transition amplitude that is associated to this algorithm,

$$\begin{array}{c}\text{ }T(ss^{})=T_0(ss^{})P_{\mathrm{acc}}(s,s^{})+\delta _{ss^{}}\delta _{ss^{}}\underset{r}{}T_0(sr)P_{\mathrm{acc}}(s,r),\text{(}\text{A}\text{.4)}\text{ }\hfill \end{array}$$

satisfies detailed balance and also conditions 1 and 3. It must be stressed, however, that for this statement to be true it is essential that $`T_0(ss^{})`$ satisfies detailed balance with respect to the distribution $`P_0(s)`$ and not only condition 2.

A.4 Using auxiliary stochastic variables

When the equilibrium probability distribution is too complicated to be treated directly, an equivalent but more accessible system may sometimes be found that involves auxiliary stochastic variables. In the case of the HMC algorithm, for example, the auxiliary variables are the pseudo-fermion field and the momenta of the link variables \[?\]. A simulation algorithm may then be set up for the enlarged system, and this leads to a correct algorithm for the original system under certain conditions.

In abstract terms the states of the enlarged system are labelled by pairs $`v,s`$, where $`v`$ and $`s`$ represent the auxiliary and the basic variables respectively. The equivalence to the original system is then expressed through the identity

$$\begin{array}{c}\text{ }P(s)=\underset{v}{}\widehat{P}(v,s)\text{(}\text{A}\text{.5)}\text{ }\hfill \end{array}$$

in which $`\widehat{P}(v,s)`$ denotes the equilibrium distribution of the enlarged system. In the following the conditional probability

$$\begin{array}{c}\text{ }\widehat{P}(v|s)=\widehat{P}(v,s)/P(s)\text{(}\text{A}\text{.6)}\text{ }\hfill \end{array}$$

to find $`v`$ given $`s`$ will play an important rôle.

Now if $`\widehat{T}(v,sv^{},s^{})`$ is any given transition probability for the enlarged system that satisfies conditions $`1`$ and $`2`$ (but not necessarily condition 3), an update algorithm for the original system is obtained as follows. Starting from an initial state $`s`$, the auxiliary variables $`v`$ are first chosen randomly with conditional probability $`\widehat{P}(v|s)`$. After that, a state $`v^{},s^{}`$ is generated with probability $`\widehat{T}(v,sv^{},s^{})`$ and $`s^{}`$ is taken to be the new state of the original system. It is straightforward to show that the associated transition probability

$$\begin{array}{c}\text{ }T(ss^{})=\underset{v,v^{}}{}\widehat{P}(v|s)\widehat{T}(v,sv^{},s^{})\text{(}\text{A}\text{.7)}\text{ }\hfill \end{array}$$

satisfies conditions 1 and 2. Moreover, detailed balance holds if the transition probability $`\widehat{T}(v,sv^{},s^{})`$ has this property relative to the equilibrium distribution of the enlarged system. Whether condition 3 is fulfilled depends on both $`\widehat{P}(v|s)`$ and $`\widehat{T}(v,sv^{},s^{})`$ and has to be verified in each case.

A.5 Stochastic acceptance--rejection method

The starting point here is again an algorithm with transition probability $`T_0(ss^{})`$ that satisfies detailed balance with respect to some approximate distribution $`P_0(s)`$. An enlarged system is then considered, exactly as above, but the auxiliary variables are now only used in the acceptance–rejection step where (A.3) is replaced by

$$\begin{array}{c}\text{ }P_{\mathrm{acc}}(s,s^{})=\underset{v}{}\widehat{P}(v|s)\mathrm{min}\{1,\frac{P_0(s)\widehat{P}(v,s^{})}{\widehat{P}(v,s)P_0(s^{})}\}.\text{(}\text{A}\text{.8)}\text{ }\hfill \end{array}$$

This means that once a new state $`s^{}`$ of the original system is proposed, a random choice of the auxiliary variables $`v`$ first needs to be made, with conditional probability $`\widehat{P}(v|s)`$, before it can be decided whether $`s^{}`$ should be accepted or not.

The correctness of this method follows from the identity

$$\begin{array}{c}\text{ }\left[P(s)/P_0(s)\right]P_{\mathrm{acc}}(s,s^{})=\left[P(s^{})/P_0(s^{})\right]P_{\mathrm{acc}}(s^{},s)\text{(}\text{A}\text{.9)}\text{ }\hfill \end{array}$$

that is easily derived from the definition (A.8) and the properties of the distributions $`P_0(s)`$ and $`\widehat{P}(v,s)`$. In particular, the associated transition probability (A.4) satisfies detailed balance. It has recently been noted, however, that \[?\]

$$\begin{array}{c}\text{ }P_{\mathrm{acc}}(s,s^{})\mathrm{min}\{\underset{v}{}\widehat{P}(v|s),\underset{v}{}\frac{P_0(s)\widehat{P}(v,s^{})}{P(s)P_0(s^{})}\}\text{ }\hfill \\ \\ \text{ }=\mathrm{min}\{1,\frac{P_0(s)P(s^{})}{P(s)P_0(s^{})}\},\text{(}\text{A}\text{.10)}\text{ }\hfill \end{array}$$

which shows that the stochastic acceptance–rejection method tends to be less efficient than the non-stochastic method.

References

 S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, Phys. Lett. B195 (1987) 216

 P. de Forcrand, T. Takaishi, Nucl. Phys. B (Proc. Suppl.) 53 (1997) 968

 R. Frezzotti, K. Jansen, Phys. Lett. B 402 (1997) 328; Nucl. Phys. B 555 (1999) 395 and 432

 M. Lüscher, Nucl. Phys. B418 (1994) 637

 H. A. Schwarz, Gesammelte Mathematische Abhandlungen, vol. 2 (Springer Verlag, Berlin, 1890)

 Y. Saad, Iterative methods for sparse linear systems (PWS Publishing Company, Boston, 1996); see also http://www-users.cs.umn.edu/~saad/

 M. Hasenbusch, Phys. Rev. D59 (1999) 054505

 R. Frezzotti et al. (ALPHA collab.), Comput. Phys. Commun. 136 (2001) 1

 A. Hasenfratz, F. Knechtli, Comput. Phys. Commun. 148 (2002) 81

 A. Alexandru, A. Hasenfratz, Phys. Rev. D65 (2002) 114506; ibid. D66 (2002) 094502

 F. Knechtli, U. Wolff, Dynamical fermions as a global correction, hep-lat/0303001

 K. G. Wilson, Phys. Rev. D10 (1974) 2445

 B. Sheikholeslami, R. Wohlert, Nucl. Phys. B 259 (1985) 572

 M. Lüscher, S. Sint, R. Sommer, P. Weisz, Nucl. Phys. B478 (1996) 365

 M. Lüscher, R. Narayanan, P. Weisz, U. Wolff, Nucl. Phys. B384 (1992) 168

 S. Sint, Nucl. Phys. B421 (1994) 135

 K. Jansen, R. Sommer (ALPHA collab.), Nucl. Phys. B530 (1998) 185 \[E: ibid. B643 (2002) 517\]

 L. Giusti, C. Hoelbling, M. Lüscher, H. Wittig, Numerical techniques for lattice QCD in the $`ϵ`$-regime, hep-lat/0212012, to appear in Comput. Phys. Commun.

 C. T. H. Davies et al., Phys. Rev. D41 (1990) 1953

 S. Aoki et al. (CP-PACS collab.), Phys. Rev. D67 (2003) 034503

 M. Hasenbusch, Phys. Lett. B519 (2001) 177

 M. Hasenbusch, K. Jansen, Speeding up lattice QCD simulations with clover improved Wilson fermions, hep-lat/0211042

 A. Duncan, E. Eichten, H. Thacker, Phys. Rev. D59 (1999) 014505

 A. Duncan, E. Eichten, Y. Yoo, Unquenched QCD with light quarks, hep-lat/0209123

 M. Peardon, J. Sexton (TrinLat Collab.), Multiple molecular dynamics time scales in hybrid Monte Carlo fermion simulations, hep-lat/0209037

 N. Metropolis et al., J. Chem. Phys. 21 (1953) 1087


format: text

text:
# The phase shift effective range expansion from supersymmetric quantum mechanics

## I Introduction

In the description of the interaction between composite particles by local potentials an ambiguity arises between different phase equivalent or isophase potential families. There are shallow potentials, which possess only physical bound states and deep potentials, which possess both physical and unphysical bound states. The latter, called Pauli forbidden states, simulate nonlocal aspects of the potential, or else, the complexity of the interaction between composite particles. The number of Pauli forbidden states can be predicted from a microscopic description of the interacting particles SAITO . The application of the supersymmetric (SUSY) quantum mechanics WITTEN  to the inverse scattering problem provides an elegant and powerful algebraic way to understand the relation between such phase equivalent potentials SUKUMAR ; BS1 ; LEEB ; SS . A supersymmetric transformation can be seen as a specific Darboux transformation. In the following we shall use Darboux and SUSY transformations as synonyms.

The above procedures cannot directly provide a correct behavior of the phase shift at energies small relative to the potential strength, i. e. a correct effective range expansion for $`\mathrm{}>0`$. In our opinion the reason is that the role of the angular momentum for a given central potential was not properly understood so far. In the framework of SUSY quantum mechanics the problem has been raised by Sukumar SUKUMAR  and Baye and Sparenberg BS1  but not solved in principle. However, it has been tackled pragmatically in Ref. SPAR . For example, in the case of the $`\mathrm{}=4`$ partial wave it led to an S-matrix containing powers of $`k`$ restricted to $`n10`$. How one arrives at such a restriction is neither explained nor the power 10 is justified.

The basic idea of this study is to associate the angular momentum with the long-distance asymptotic behavior of the potential, irrespective of its singularity at the origin. This is in the spirit of Ref. SWAN  where these asymptotic limits are independent of each other. This starting point will provide a new possibility for getting a correct effective range expansion of the phase shift which is the following Taylor series expansion in the vicinity of $`k=0`$ (see e. g. NEWTON )

$$k^{2\mathrm{}+1}\mathrm{cot}\delta _{\mathrm{}}(k)=\frac{1}{a_0\mathrm{}}+\frac{1}{2}r_0\mathrm{}k^2+\mathrm{}.$$
(1)
Here $`a_0\mathrm{}`$ is the scattering length and $`r_0\mathrm{}`$ the effective range. The expression (1) implies that for a given $`\mathrm{}`$, in the series expansion of $`\mathrm{tan}\delta _{\mathrm{}}(k)`$ the coefficients of the terms containing powers of $`k`$ below $`2\mathrm{}+1`$ must vanish. In the frame of SUSY quantum mechanics we solve this problem by introducing adequate long-distance Darboux transformations.

The paper is organized as follows. In the subsection A of the next section we introduce the Darboux transformation method and briefly review the $`\mathrm{}`$-fixed transformations. In subsection B we introduce the long-distance transformations. Section III is devoted to results and applications to the neutron-proton scattering. Details are worked out for the $`\mathrm{}`$=1 and $`\mathrm{}`$=2 partial waves. Conclusions are drawn in the last section.

## II Theory

### II.1 $`\mathrm{}`$-fixed Darboux transformations

We recall that the Darboux transformation method consists in getting solutions $`\phi `$ of one Schrödinger equation

$$h_1\phi =E\phi ,h_1=\frac{d^2}{dx^2}+V_1(x),$$
(2)
when solutions $`\psi `$ of another equation

$$h_0\psi =E\psi ,h_0=\frac{d^2}{dx^2}+V_0(x),$$
(3)
are known. This is achieved by acting on $`\psi `$ with a differential operator $`L`$ of the form

$$\phi =L\psi ,L=d/dx+w(x),$$
(4)
where the real function $`w(x)`$, called superpotential, is defined as the logarithmic derivative of a known solution of (3) denoted by $`u`$ in the following. One has

$$w=u^{}(x)/u(x),h_0u=\alpha u,$$
(5)
with $`\alpha E_0`$, where $`E_0`$ is the ground state energy of $`h_0`$ if it has a discrete spectrum or the lower bound of the continuous spectrum otherwise. The function $`u`$ is called transformation or factorization function and $`\alpha `$ its factorization constant or factorization energy. The potential $`V_1`$ is defined in terms of the superpotential $`w`$ as

$$V_1(x)=V_0(x)2w^{}(x).$$
(6)
Eq. (4) defines a first order Darboux transformation. In the following we shall deal with chains of $`N`$ successive transformations of this type.

Let us start by first considering $`\mathrm{}`$-fixed transformations as in Ref. SS . This means that we use a special chain of $`N=2n\nu `$ first order Darboux transformations with $`\nu 0`$, generated by the following system of transformation functions

$`v_1(x),\mathrm{},v_\nu (x),u_{\nu +1}(x),v_{\nu +1}(x),\mathrm{},u_n(x),v_n(x)`$ (7)
$`h_0u_j(x)=a_j^2u_j(x),h_0v_j(x)=b_j^2v_j(x),`$ (8)
where $`v_j`$ are regular ($`v_j(0)=0`$) and $`u_j`$ irregular ($`u_j(0)0)`$, the latter being expressed in terms of the Jost solutions as

$$u_j(x)=A_jf(x,ia_j)+B_jf(x,ia_j).$$
(9)
They have arbitrary eigenvalues $`a_{j}^{}{}_{}{}^{2}`$ and $`b_{j}^{}{}_{}{}^{2}`$ respectively, but always below $`E_0`$. If we are interested in the final action of the chain only, the solution $`\psi _N(x,k)`$ of the transformed equation with the Hamiltonian

$$h_N=d^2/dx^2+V_N$$
(10)
corresponding to the energy $`E=k^2`$ is given by CRUM

$$\psi _N(x,k)=W(u_1,\mathrm{},u_N,\psi _0(x,k))W^1(u_1,\mathrm{},u_N)$$
(11)
where $`W`$ are Wronskians expressed in terms of $`u_j`$, denoting symbolically any function of (7) and of $`\psi _0(x,k)`$ which is a solution of the original Schrödinger equation corresponding to the same energy $`E`$. In the Hamiltonian (10) the transformed potential is

$$V_N=V_02\frac{d^2}{dx^2}\mathrm{ln}W(u_1,\mathrm{},u_N).$$
(12)
For $`N=1`$ one has $`W(u_1)u_1`$ and one recovers (6) with $`u=u_1`$. If $`V_0`$ is finite at the origin, $`V_N`$ behaves as $`\nu (\nu +1)x^2`$ when $`x0`$. Therefore the parameter $`\nu `$ is called the singularity strength. The formulas (11) and (12) result from the replacement of a chain of $`N`$ first order transformations by a single $`N`$th order transformation, which happens to be more efficient in practical calculations.

In Ref. SS  we obtained that the transformed Jost function $`F_N`$ is related to the initial Jost function $`F_0`$ by

$$F_N(k)=F_0(k)\underset{j=\mathrm{\hspace{0.17em}1}}{\overset{\nu }{}}\frac{k}{k+ib_j}\underset{j=\nu +1}{\overset{n}{}}\frac{kia_j}{k+ib_j}.$$
(13)
For $`\nu =0`$ the first product is unity. Since a Jost function is analytic in the upper half of the complex $`k`$-plane (see e. g. FADDEEV ) all $`b`$’s must be positive whereas the $`a`$’s can have any sign, so that every positive $`a_j`$ corresponds to a discrete level $`E=a_j^2`$ of $`h_N`$.

The corresponding phase shift $`\delta _{\mathrm{}}^N(k)`$ can be written as

$$\delta _{\mathrm{}}^N(k)=\delta _{\mathrm{}}^{\mathrm{\hspace{0.17em}0}}(k)+\mathrm{\Delta }_{\mathrm{}}^N(k),$$
(14)
where $`\delta _{\mathrm{}}^0(k)`$ is the initial phase shift due to the potential $`V_0`$ and $`\mathrm{\Delta }_{\mathrm{}}^N(k)`$ is the phase shift produced by the chain of $`N`$ Darboux transformations

$$\mathrm{\Delta }_{\mathrm{}}^N(k)=\underset{j=\nu +1}{\overset{n}{}}\mathrm{arctan}(k/a_j)\underset{j=\mathrm{\hspace{0.17em}1}}{\overset{n}{}}\mathrm{arctan}(k/b_j).$$
(15)
This is consistent with the asymptotic form of the scattering solution $`\mathrm{sin}(kx\frac{\pi \mathrm{}}{2}+\delta _{\mathrm{}}^N)`$. In the limit $`k\mathrm{}`$ one has $`\delta _{\mathrm{}}^N(\mathrm{}\nu )\pi /2`$, in agreement with Ref. SWAN  for a singular potential of parameter $`\nu `$. More detailed discussion of properties of $`\mathrm{}`$-fixed transformations may be found in SS .

In the case $`\mathrm{}=0`$ by expanding $`\delta _{\mathrm{}}^0(k)`$ and the arc tangent functions in (15) in power series one obtains the effective range expansion (1). For $`\mathrm{}>0`$ the situation is more subtle, since the first term in power series of arc tangent functions is proportional to $`k`$. Therefore in the usual practice based on the SUSY approach, where $`\delta _{\mathrm{}}^0(k)`$ is fixed, it would be difficult to cancel the undesired powers of $`k`$ in order to comply with (1). As mentioned above, we believe that the reason is that one deals with Darboux transformations which do not affect the long distance behavior of the resulting potential, as it was pointed out in Ref. SS .

### II.2 Long-distance Darboux transformations

To change the long distance behavior of a potential by SUSY transformations we use transformation functions with zero eigenvalue. To cancel undesired powers in the series expansion of arc tangent functions we derive a proper background phase shift as shown below.

Consider the potential $`V_0(x)`$ which for $`x\mathrm{}`$ behaves as

$$V_0(x)=\frac{\mathrm{}(\mathrm{}+1)}{x^2}+O\left(x^{2\gamma }\right),\gamma >0,\mathrm{}0.$$
(16)
As it is known FADDEEV  the Schrödinger equation containing a potential satisfying (16) has zero eigenvalue solutions with the following asymptotic behavior at $`x\mathrm{}`$

$`v(x)`$ $`=`$ $`Cx^{\mathrm{}+1}\left[1+O\left(x^\gamma \right)\right]`$ (17)
$`u(x)`$ $`=`$ $`{\displaystyle \frac{D}{x^{\mathrm{}}}}\left[\mathrm{\hspace{0.17em}1}+O\left(x^\gamma \right)\right].`$ (18)
The functions (17) are regular at the origin but singular at infinity and the functions (18) are just the other way round. When these functions are taken as transformation functions the change in the potential for sufficiently large $`x`$ has one of the following forms

$`\mathrm{\Delta }V(x)`$ $`=`$ $`2[\mathrm{ln}v(x)]^{\prime \prime }={\displaystyle \frac{2(\mathrm{}+1)}{x^2}}+o(x^{2\gamma })`$ (19)
$`\mathrm{\Delta }V(x)`$ $`=`$ $`2[\mathrm{ln}u(x)]^{\prime \prime }={\displaystyle \frac{2\mathrm{}}{x^2}}+o(x^{2\gamma })`$ (20)
where $`\mathrm{\Delta }V(x)=V_NV_0`$ with $`N`$ = 1 in (12). It is clear from here that the function (17) increases the value of $`\mathrm{}`$ by one unit and the function (18) decreases it by one unit. Moreover a linear combination of (17) and (18) is a function of type (17). They form a one-parameter family, while the function (18) is uniquely defined (up to an inessential constant factor). In this family there is only one function regular at the origin. This function, used as transformation function in the Darboux algorithm, changes both $`\nu `$ and $`\mathrm{}`$ but all the other members of the singular family change $`\mathrm{}`$ without affecting $`\nu `$. In the following we shall use the singular functions of the one-parameter family defined above to derive phase shifts leading to a correct effective range expansion. We shall therefore show that the parameters appearing in the linear combination of (17) and (18) can be chosen such as the resulting phase shift provides the general effective range expansion (1). Hence, starting with a given $`V_0`$ with $`\mathrm{}=0`$ we first perform a number $`N=\mathrm{}`$ of transformations which give the correct long distance behavior of the potential $`V_{\mathrm{}}`$ and introduce $`\mathrm{}`$ parameters in the phase shift. Next, an $`\mathrm{}`$-fixed chain is performed, producing the final phase shift (14) for which the potential $`V_{\mathrm{}}`$ plays the role of the initial potential. The latter transformation does not affect the asymptotic form of the potential $`V_{\mathrm{}}`$ at large $`x`$. Hence, the resulting potential $`V_N`$ has an asymptotic behavior corresponding to the $`\mathrm{}`$th partial wave. The addition of $`\mathrm{}`$ zero-energy eigenfunctions to the $`N=2n\nu `$ first order transformation functions used in the $`\mathrm{}`$-fixed chain increases $`N`$ by $`\mathrm{}`$ units. This means that in the $`N`$th order transformation to be used below the total number of transformation functions is

$$N=2n\nu +\mathrm{}.$$
(21)
If we start with the zero initial potential, $`V_00`$, the formula (14) for the phase shift has to be modified as follows:

$$\delta _{\mathrm{}}^N(k)=\delta _{\mathrm{}}^{\mathrm{}}(k)+\mathrm{\Delta }_{\mathrm{}}^N\mathrm{}(k).$$
(22)
Here $`\delta _{\mathrm{}}^{\mathrm{}}(k)`$ is produced by the long-distance transformations which give rise to an intermediate potential $`V_{\mathrm{}}`$. In the following $`\delta _{\mathrm{}}^{\mathrm{}}(k)`$ will play the role of a background phase shift. The additional phase shift $`\mathrm{\Delta }_{\mathrm{}}^N\mathrm{}(k)`$ corresponding to the $`\mathrm{}`$-fixed subchain of $`N\mathrm{}=2n\nu `$ transformations has the same form as (15). We shall illustrate this procedure by applications given in the following section.

## III Applications

### III.1 The case $`\mathrm{}=1`$

We start with the potential $`V_0=0`$ in Eq. (3). Let us take $`u_{0,1}=x+x_0`$ as the transformation function with zero eigenvalue, where $`x_00`$ is a free parameter. Then the first order transformation operator (4) takes the form

$$L_1=\frac{d}{dx}+\frac{1}{x+x_0}.$$
(23)
The transformed potential is

$$V_1=\frac{2}{(x+x_0)^2},x_00$$
(24)
and its Jost solution $`f_1(x,k)`$ may be found by applying the operator (23) on the Jost solution $`f_0(x,k)=\mathrm{exp}(ikx)`$ of the free particle equation. After dividing by the factor $`ik`$ one finds

$$f_1(x,k)=\left(1\frac{1}{ik(x+x_0)}\right)\mathrm{exp}(ikx).$$
(25)
For $`x_00`$ the potential (24) has $`\nu =\nu _1=0`$ and its Jost function $`F_1(k)`$ coincides with $`f_1(0,k)`$. If one now applies the operator (23) on an oscillating solution of the free particle equation $`\mathrm{sin}(kx+\delta _1^1)`$ one obtains

$$\phi _1(x,k)=k\mathrm{cos}(kx+\delta _1^1)+\frac{\mathrm{sin}(kx+\delta _1^1)}{x+x_0}.$$
(26)
This solution is regular at the origin provided

$$\delta _1^1=\mathrm{arctan}kx_0,$$
(27)
and has the asymptotic behavior $`\mathrm{sin}(kx+\delta _1^1\frac{\pi }{2})`$ at $`x\mathrm{}`$, it describes the $`\mathrm{}=1`$ scattering state of the potential (24). Thus we have switched from the partial wave $`\mathrm{}=0`$ of $`V_0=0`$ to $`\mathrm{}=1`$ of the potential (24). Now we can perform $`\mathrm{}`$-fixed transformations. The free parameter $`x_0`$ will be chosen so that the final phase shift (22) will have the correct effective range expansion. Replacing $`\delta _{\mathrm{}}^{\mathrm{}}`$ by the result (27) and expanding all arc tangent functions in power series one can see that the coefficient of the term proportional to $`k`$ in $`\mathrm{tan}\delta _1^N`$ vanishes for

$$x_0=\underset{j=\nu +1}{\overset{n}{}}a_j^1+\underset{j=1}{\overset{n}{}}b_j^1.$$
(28)
Both the regular and irregular solutions corresponding to the potential (24) can be found with the help of the Jost solution (25). But in the spirit of Ref. CRUM , we can avoid this step, thus considerably reducing the amount of numerical work. This means that in the formulas (11) and (12) we can directly use appropriate solutions of the free particle equation which are simple linear combinations of exponentials

$$\psi _j(x,k)=A_j\mathrm{exp}(ik_jx)+B_j\mathrm{exp}(ik_jx),\text{Re}k_j=0,$$

where we have to find the correct ratio $`B_j/A_j`$. Since the regular solutions of the potential (24), consisting of functions $`v_j(x)=L_1\psi _j(x,ib_j)`$ satisfy the condition $`v_j(0)=0`$, which fixes the ratio $`B_j/A_j`$, we have free particle solutions of the form

$$\psi _j(x,ib_j)=(b_jx_0+1)\mathrm{exp}(b_jx)+(b_jx_01)\mathrm{exp}(b_jx).$$
(29)
The irregular solutions for the same potential, defined as $`u_j=L_1\psi _j(x,ia_j)`$, should be obtained from the functions

$$\psi _j(x,ia_j)=A_j\mathrm{exp}(a_jx)+B_j\mathrm{exp}(a_jx),$$
(30)
which for $`B_j/A_j(a_jx_01)/(a_jx_0+1)`$, $`A_j0`$ and $`a_j>0`$ have an increasing asymptotic behavior but if $`A_j=0`$ they decrease asymptotically. To stress the difference between solutions with different asymptotic behavior we choose in the latter case $`a_j<0`$ and $`A_j0`$ but $`B_j=0`$ (for more details see SS ). Here $`a_j`$, $`b_j`$ and $`x_0`$ are parameters of the model to be found below. Then in the $`N`$th order transformation we use the functions $`u_{0,1}=x+x_0`$, (29) and (30) to calculate (12). Note nevertheless that in the resulting phase shift, given by (22), $`\delta _{\mathrm{}}^{\mathrm{}}`$ has to be replaced by $`\delta _1^1`$ of (27) since the initial potential for the subchain of $`\mathrm{}`$-fixed transformations is now $`V_1`$ of (24). For the $`N`$th order transformation the potential $`V_1`$ and the phase shift $`\delta _1^1`$ play an auxiliary role.

As an application we look for neutron-proton ($`np`$) potentials which reproduce the ’pruned’ phase shift of STOKS  for the $`{}_{}{}^{1}P_{1}^{}`$ partial wave. This phase shift together with the theoretical values obtained from the expression (14) and denoted by $`\delta _1^7`$ are exhibited in Table 1. The six fitted $`S`$-matrix poles are

$$\begin{array}{c}a_1=0.7290,a_2=0.7295,a_3=1.0368\hfill \\ b_1=0.4403,b_2=0.4408,b_3=3.3818\hfill \end{array}$$
(31)
in fm<sup>-1</sup> units. The superscript 7 carried by the phase shift is consistent with the formulas (21) and (22) and it implies that we used a $`7`$th order transformation according to (11) and (12). Then from Eq. (28) one gets $`x_0=3.0578`$ fm. Now we can expand all seven arc tangent functions appearing in (14) in power series. This leads to a correct effective range expansion given by

$$k^3\mathrm{cot}\delta _1^7(k)=0.31823.1511k^2+\mathrm{}$$
(32)
from which one can extract the scattering length $`a_{01}`$ and the effective range $`r_{01}`$ defined according to (1). In Table 2 these values are compared with another theoretical model RIJKEN . They are surprisingly close to each other. The phenomenological Reid93 potential STOKS  also gives similar values MART . Moreover the scattering length is located in the interval deduced from a partial wave analysis MATHELITSCH .

In order to construct potentials giving rise to the phase shift $`\delta _1^7`$, the poles (31) have to be associated with transformation functions defined by (29) and (30). The poles $`b_j`$ correspond to regular solutions of $`V_1(x)`$ resulting from (29). The poles $`a_1`$ and $`a_2`$, which are negative, correspond to decreasing functions of the form (30) with $`B_j=0`$ ($`j`$=1,2). It remains the pole $`a_3`$ which is positive. If we take $`B_3/A_3(a_3x_01)/(a_3x_0+1)`$ in Eq. (30), we obtain from (12) a one-parameter ($`B_3/A_3`$) family of one-level isophase deep potentials with the discrete level at $`E=a_3^2`$. But if we choose $`B_3/A_3=(a_3x_01)/(a_3x_0+1)`$ the initially irregular function moves into the regular family,

the level disappears and one gets the uniquely defined shallow representative of the family of isophase potentials, denoted by $`V_7`$. Fig. 1 shows this potential from which the centrifugal term has been subtracted. The potential $`V_7`$, is quite close to the Reid68 potential REID68 , represented in the same figure. Figure 2 shows one of the deep isophase potentials corresponding to $`B_3/A_3=0.351`$. This constant has been adjusted to get a potential as close as possible to Reid93 STOKS  in the interval $`0.4`$ fm $`x2`$ fm. This deep potential possesses a Pauli forbidden state of energy $`E=a_3^2=44.46`$ MeV. Contrary to the Reid93 potential which is deep but finite, our potential behaves at origin as $`2/x^2`$. If necessary, it can be regularized as for example in Ref. SPAR .

It would be interesting to analyze the asymptotic behavior of our calculated potential to see if it is compatible with modern phenomenological potentials constructed in the spirit of meson theory, i.e. which include one-pion-exchange (OPE) contributions. This is precisely the case of the Reid93 potential, which includes OPE with neutral and charged pions (for details see Ref. STOKS ). The comparison given in Figure 3 shows that the shallow potential $`V_7`$ and its deep partner are practically identical in the asymptotic region and that they are both extremely close to the asymptotic form of the updated Reid93 potential MART . Thus our potentials are in the excellent agreement with expectations from Yukawa’s OPE theory. We remind that we subtracted the centrifugal barrier from our potentials in order to make this comparison feasible. In solving the Schrödinger equation this should be added back to the nucleon-nucleon interaction in each case. The fact that the shallow and deep potentials are asymptotically the same is physically correct. As mentioned in the introduction, deep potentials reflect the compositeness (quark structure) of the interacting particles in the overlap region but, once the particles are well separated, they could be treated as point particles as in OPE theories, which means that beyond some distance the deep and shallow potentials should coincide.

### III.2 The case $`\mathrm{}=2`$.

Now we have to apply two subsequent transformations associated with functions corresponding to $`E=0`$ and then, as above, a subchain of $`\mathrm{}`$-fixed transformations. After the first transformation with the function $`u_{0,1}=x+x_0`$, the potential $`V_1`$ of (24) has $`u=1/(x+x_0)`$ and $`\stackrel{~}{u}=(x+x_0)^2`$ as linearly independent solutions at $`E=0`$. Their linear combination $`u_{0,2}=cu+\stackrel{~}{u}`$, which is the transformation function for the second transformation step defined by the operator

$$L_2=d/dx+u_{0,2}^{}(x)/u_{0,2}(x)$$
(33)
contains two free parameters $`x_0`$ and $`c`$. These can be chosen such as the series expansion for $`\mathrm{tan}\delta _2(k)`$ starts at $`k^5`$. The intermediate (or background) potential

$$V_2=2[\mathrm{ln}u_{0,2}(x)u_{0,1}(x)]^{\prime \prime }=\frac{6(x+x_0)[(x+x_0)^32c]}{[c+(x+x_0)^3]^2}$$
(34)
obtained from the $`L=L_2L_1`$ Darboux transformation (for more details see SS ) plays now the role of the initial potential for an $`\mathrm{}`$-fixed subchain of transformations. The background phase shift (modulo $`\pi `$) corresponding to $`V_2`$ is

$$\delta _2^2=\mathrm{arctan}\frac{3kx_0^2}{3x_0k^2(x_0^3+c)}.$$
(35)
Note that the function $`\phi _2(x,k)=L_2L_1\mathrm{sin}(kx+\delta _2^2)`$ is regular at the origin and describes an unnormalized scattering state for $`\mathrm{}=2`$. In the formula (22) we have to identify $`\delta _{\mathrm{}}^{\mathrm{}}(k)`$ with $`\delta _2^2(k)`$ of (35). After expanding in power series all arc tangent functions we find that the coefficient of the term linear in $`k`$ vanishes for $`x_0`$ given by (28) and

$$c=\frac{1}{3}[\underset{j=1}{\overset{n}{}}(a_j^3+b_j^3)]$$
(36)
ensures the cancellation of the coefficient of $`k^3`$ in the series expansion of $`\mathrm{tan}\delta _2^N`$. Now, to find solutions of the free particle equations, giving rise to the regular family of the potential (34), to be used in (12), we need eigenfunctions of $`h_0=d^2/dx^2`$ satisfying the condition $`L_2L_1\psi _j(x,ib_j)=0`$ at $`x=0`$. They are given by the following linear combination of exponentials

$`\psi _j(x,ib_j)`$ $`=`$ $`[\mathrm{\hspace{0.17em}3}x_0+3b_jx_0^2+b_j^2(c+x_0^3)]\mathrm{exp}(b_jx)`$
$``$ $`[\mathrm{\hspace{0.17em}3}x_03b_jx_0^2+b_j^2(c+x_0^3)]\mathrm{exp}(b_jx).`$
The irregular family still results from (30) subject to the condition that the ratio $`B_j/A_j`$ is different from that presented in (III.2).

With a fit of a similar quality to that performed for $`{}_{}{}^{1}P_{1}^{}`$ we could reproduce the $`{}_{}{}^{1}D_{2}^{}`$ partial wave phase shift of STOKS  with the following four poles of the $`S`$-matrix

$$\begin{array}{c}a_1=0.2047,a_2=1.9800\hfill \\ b_1=1.2305,b_2=4.9631\hfill \end{array}$$
(38)
in fm<sup>-1</sup> units. From Eqs. (28) and (36) we get $`x_0=4.375`$ fm, $`c=116.08`$ fm<sup>3</sup>. This leads to the following effective range expansion

$$k^5\mathrm{cot}\delta _2^6(k)=0.4496+10.9878k^2+\mathrm{}.$$
(39)
where the superscript $`N=6`$ represents 4 transformation functions associated with the poles (38) plus two zero-eigenvalue functions $`u_{0,1}`$ and $`u_{0,2}`$ defined above. This is consistent with the formula (21) with $`n=2`$, $`\nu =0`$ and $`\mathrm{}=2`$. From the expansion (39) one can extract the scattering length $`a_{02}`$ and the effective range $`r_{02}`$ defined in Eq. (1). These values are shown in Table 2. They are comparable to those obtained for the potential Reid93 MART .

## IV Conclusions

By working out these two particular cases we have shown that a new insight emerges into the role of the angular momentum of a central potential. If associated with the long-distance behavior of the potential, it allows us to introduce transformations that bring $`\mathrm{}`$ free parameters in the background phase shift (Eqs (27) and (35)). When a correct effective range expansion is required, each extra unit of angular momentum imposes a new constraint on the whole system of parameters of the model such as the number of constraints coincides with the number of parameters in the background phase shift. For the particular cases of $`\mathrm{}=1`$ and $`\mathrm{}=2`$ explicit solutions of the constraint equations are given. Thus, the extension of the method to $`\mathrm{}>2`$ is straightforward.

Before ending we should mention that the generalized Levinson theorem SWAN  is always satisfied in our approach.

We are grateful to Mart Rentemeester for useful correspondence. B.F.S. acknowledges support from the Spanish Ministerio de Education, Cultura y Deporte Grant SAB2000-0240 and the Spanish MCYT and European FEDER grant BFM2002-03773. He is also grateful for hospitality at the Fundamental Theoretical Physics Laboratory of the University of Liege.