text:
# Parallel algorithm with spectral convergence for nonlinear integro-differential equations

## 1 Introduction

Astrophysical applications related to the physics of the early universe, as well as challenges posed by the physics programs at new heavy ion accelerators, have triggered a renewed interest in the understanding of real time processes in the context of quantum field theory. With the advent of new computer technology and the recent success of new computational schemes, non-equilibrium phenomena which have been previously studied only in the framework mean-field theory , are now being revisited, and more complex next to leading order approaches  are being used in an attempt to clarify the role played by the rescattering mechanism, which is responsible for driving an out of equilibrium system back to equilibrium. Of particular interest is the study of the dynamics of phase transitions and particle production following a relativistic heavy-ion collision. One way of approaching this study is based on solving Schwinger Dyson equations within the closed time path (CTP) formulation . This formalism has been recently shown to provide good approximations of the real time evolution of the system both in quantum mechanics and 1+1 dimensional classical field theory , where direct comparisons with exact calculations can be performed.

The key element in carrying out such studies is related to the calculation of the two-point Green function, which is solved for self-consistently with the equations of motion for the fields. The two-point Green function gives rise to Volterra-type integral or integro-differential equations. In the process of extending our study to encompass a higher number of spatial dimensions, i.e. 2+1 and 3+1 field theory, we are faced with the challenge of coping with constraints dictated both by storage and time-related computational limits. Thus our interest in designing algorithms which feature spectral convergence in order to achieve convergence with minimum storage requirements. In addition, we also desire these algorithms to scale when ported to massively multiprocessor (MPP) machines, so that solutions can be obtained in a reasonable amount of time.

Algorithms for Volterra integral and integro-differential equations usually start out at the lower end of the domain, $`a`$, and march out from $`x=a`$, building up the solution as they go . Such methods are *serial* by nature, and are, in general, not suitable for parallel implementation on a MPP machine. Even so, clever approaches to already existing methods can provide algorithms that take advantage of a parallel processing computer: Shaw  has shown recently that once the starting values of the approximation are obtained, one can design a *global* approach where successive approximations of the solution over the entire domain $`x[a,b]`$ can be evaluated simultaneously.

In a recent paper  one of us has discussed a spectral method  of solving some types of equations of interest for the study of time-dependent nonequilibrium problems in quantum field theory. The gist of the method consists in expanding out the unknown function in terms of Chebyshev polynomials on a suitable grid, thus reducing the problem to finding the numerical solution of a system of linear equations. The main advantage of this method over standard finite-difference type methods resides in the spectral character of its convergence. This is related in part to the fact that Chebyshev type methods use a non-uniform grid, while finite-difference methods require a uniform grid. Usually there is a trade-off between computational time and storage requirements, and a balanced solution must be reached on a case-by-case basis. Spectral methods are more expensive per point as the matrices may be considerably denser than in the finite-difference case, but we require considerably fewer grid points in order to achieve the same degree of accuracy. By expanding the unknown function on a compact support in Chebyshev polynomials and using a partition of the domain based either on the set of (N+1) extrema or the set of N zeros of $`T_N(x)`$ – the Chebyshev polynomial of first kind of degree N – we in fact replace a continuous problem by a discrete one. For non-singular functions the discrete orthogonality and completeness relations for Chebyshev polynomials at the above grid points assure a *defacto* exact expansion for an *arbitrary* finite value $`N`$. In practice however, one has to compute derivatives and integrals of the unknown function at the collocation points, and the Chebyshev expansion provides only an approximation for these subsequent computations. These errors, together with the finite accuracy of numerical methods needed in conjunction with the Chebyshev expansion, conspire in order to deteriorate the accuracy of the solution at very small values of N.

The paper is organized as follows: In Section 2, for comparison purposes, we start by reviewing a finite-difference approach for the numerical solution of Volterra type integro-differential equations. We review the general framework of the Chebyshev-expansion method in Section 3, and illustrate our approach for the case of Volterra integro-differential equations. In Section 4 we present a complete assessment of the convergence and computational cost of the proposed method for the case of a test problem, and compare with results obtained via the finite-difference method. In Section 5 we discuss the relevant aspects of a large-scale calculation arising in the study of time-dependent quantum field theory, for which our numerical strategy is particularly suitable. We present our conclusions in Section 6.

## 2 Stable multi-step method for Volterra type equations

The type of problems arising in the study of time-dependent nonequilibrium quantum field theory via a Schwinger-Dyson equation approach, can be formally reduced to the general case of a nonlinear Volterra integro-differential equation. Direct methods for solving nonlinear Volterra integral and integro-differential equations are inherently serial and therefore have not received much attention for use on a parallel computer. It is worth mentioning here the work of Crisci *et al* , who concentrated on the stability aspects of parallel iteration of Volterra-Runge-Kutta (VRK) methods for solving Volterra integral equations on parallel computers. VRK methods are step-by-step methods and can take advantage of parallel architecture. Sommeijer *et al*  covered the stability of parallel block methods for ordinary differential equations (ODE) and included equations of the integro-differential type in their discussion.

We summarize here a recent parallel algorithm , which concentrates on modifying the algorithmic side of the numerical solution process for use on a parallel processor while consciously utilizing methods that are known to be stable. The algorithm is in effect an example of a higher-order finite-difference approach, and we use this approach to compare with the spectral method presented later in this paper.

For illustration, let us consider a first-order nonlinear Volterra integro-differential equation of the form

$`𝐲^{}(x)=F[x,𝐲,Z[x;𝐲]],x[a,b],`$ (1)
with

$`Z[x;𝐲]={\displaystyle _a^x}K[x,t;𝐲(t)]dt,`$ (2)
and subject to the initial condition

$`𝐲(a)=y_0.`$ (3)
Let $`I_N`$ be a partition of I=\[a,b\], where $`I_N=\{x_N=a+nh,n=0(1)N,Nh=(ba)\}`$. The problem is to find approximations $`y_n`$ to the solution $`𝐲(x_n)`$ of Eqs. (13) for each $`x_nI_N`$. A $`k`$-step method for an integro-differential equation of the form (1) is given by

$`y_{n+1}=y_n+h{\displaystyle \underset{j=0}{\overset{k}{}}}w_jF(x_{nj},y_{nj},z_{nj}),n=k(1)N,`$ (4)
where

$`z_{nj}=h{\displaystyle \underset{i=0}{\overset{nj}{}}}c_{nj,i}K(x_{nj},x_i,y_i),j=0(1)k,y_0=𝐲(a).`$ (5)
The weights $`w_i`$ depend on the $`k`$-step method selected and the weights $`c_{i,j}`$ are those of a standard quadrature formula for integrating a function whose value is known at equally spaced steps, such as a Newton-Cotes or Newton-Gregory quadrature rule. For our multi-step ($`k=4`$) method  we choose the fourth order Adams-Bashforth predictor

$`y_{k+1}^0`$ $`=y_k+{\displaystyle \frac{h}{24}}[`$ $`55F(x_k,y_k,z_k)59F(x_{k1},y_{k1},z_{k1})`$
$`+37F(x_{k2},y_{k2},z_{k2})9F(x_{k3},y_{k3},z_{k3})],`$
and the Adams-Moulton corrector

$`y_{k+1}`$ $`=y_k+{\displaystyle \frac{h}{24}}[`$ $`9F(x_{k+1},y_{k+1}^0,z_{k+1})+19F(x_k,y_k,z_k)`$
$`5F(x_{k1},y_{k1},z_{k1})+F(x_{k2},y_{k2},z_{k2})],`$
while the integral term (2) is calculated based on the Newton-Gregory quadrature formula. We use a fourth order Runge-Kutta method in order to start out the calculation.

In order to make the algorithm suitable for parallel processing, it is useful to recall that a standard quadrature method based on an uniform grid for the integral term $`z_i`$ requires knowledge of the integrand function at the abscissas in the interval $`[x_0,x_i]`$. This is obviously a serial process and not a good candidate for parallelization. It can be observed however, that once the starting values are obtained, *all* approximations $`z_i`$ with $`i=0(1)k1`$ can simultaneously be evaluated up to and including $`x_{k1}`$. After that, once a value of $`y_j`$ corresponding to a new step $`x_j`$ is established via the predictor-corrector method, all values $`z_i`$ with $`i=j(1)N`$ can also be evaluated simultaneously. This observation makes the following algorithm possible:

1. Find the starting values $`(y_i,z_i)`$ with $`i=0(1)k1`$
2. $`doi=k,N`$

add contributions to $`z_i`$ corresponding to $`(x_j,y_j)`$, where $`j=0(1)k1`$
3. $`doi=k,N`$

1. predict $`y_i`$
    2. estimate $`z_i`$ from $`(x_i,y_i)`$
    3. correct $`y_i`$
    4. $`doj=i,N`$

    update $`z_j`$ by adding the contribution corresponding to $`(x_i,y_i)`$

The above numerical algorithm is implemented using the OpenMP style directives for the Portland Group’s pgf77 FORTRAN compiler, and reportedly shows good scalability on a shared-memory multiprocessor. The speedup of the finite difference method is best for a large number of grid points which, correspondingly, gives a better solution approximation. For example, with N=5120 and 4 processors the speedup is 3.86,a good measure of processor utilization.

While the preceding algorithm performs well on a shared memory platform, it does not port easily to an MPP machine. Before we comment on the efficiency of the algorithm, let us make two general comments: Firstly, we denote by $`T_{\text{calc}}`$ and $`T_{\text{comm}}`$ the time required to perform a floating-point operation and the time required to send a floating-point number, respectively. Secondly, we will ignore for simplicity the effect of message sizes on communication costs, and assume throughout that the ratio $`T_{\text{comm}}/T_{\text{calc}}`$ is independent of $`N`$.

Returning now, to our proposed algorithm, we remark that the communication cost for the corresponding implementation involves only the integral terms. Even so, using the message-passing interface (MPI) protocol the communication cost is $`4\mathrm{log}N`$ for the starting values and up to $`N^2`$ for the remainder of the algorithm which gives a total of $`(N^2+4\mathrm{log}N)T_{\text{comm}}`$. The total number of flops depends on the specific application but a reasonable measure is the number of function evaluations which is given by $`(N^2+4N)T_{\text{calc}}`$. The ratio of communication to computation

$$\frac{N^2+4\mathrm{log}N}{N^2+4N}\frac{T_{\text{comm}}}{T_{\text{calc}}}$$

approaches a *constant* value as $`N`$ gets larger. The communication overhead problem can be relaxed by employing a spectral method discussed in the following section, the improvement being especially significant for a multidimensional problem of the type required by our nonequilibrium quantum field theory calculations .

## 3 Spectral method with Chebyshev polynomials

Consider the $`N+1`$ extrema of the Chebyshev polynomial of the first kind of degree $`N`$, $`T_N(x)`$. This set defines a non-uniform grid in the interval $`[1,1]`$, as

$$\stackrel{~}{x}_k=\mathrm{cos}\left(\frac{\pi k}{n}\right),k=0(1)N.$$
(8)
On this grid, the Chebyshev polynomials of degree $`i<n`$ obey discrete orthogonality relations

$$\underset{k=0}{\overset{N}{}}{}_{}{}^{\prime \prime }T_{i}^{}(\stackrel{~}{x}_k)T_j(\stackrel{~}{x}_k)=\beta _i\delta _{ij},$$
(9)
where the constants $`\beta _i`$ are

$$\beta _i=\{\begin{array}{cc}\frac{N}{2},\hfill & i0,N,\hfill \\ N,\hfill & i=0,N.\hfill \end{array}$$

Here, the summation symbol with double primes denotes a sum with both the first and last terms halved. We approximate an arbitrary continuous function of bounded variation $`f(x)`$ in the interval $`[1,1]`$, as

$`f(x)`$ $``$ $`{\displaystyle \underset{j=0}{\overset{N}{}}}{}_{}{}^{\prime \prime }b_{j}^{}T_j(x),`$ (10)
with

$`b_j`$ $`=`$ $`{\displaystyle \frac{2}{N}}{\displaystyle \underset{k=0}{\overset{N}{}}}{}_{}{}^{\prime \prime }f(\stackrel{~}{x}_k)T_j(\stackrel{~}{x}_k),j=0(1)N.`$ (11)
Eq. (10) is exact at x equal to $`\stackrel{~}{x}_k`$ given by Eq. (8). Based on Eq. (10), we can also approximate derivatives and integrals as

$`f^{}(x)`$ $``$ $`{\displaystyle \underset{k=0}{\overset{N}{}}}{}_{}{}^{\prime \prime }f(\stackrel{~}{x}_k){\displaystyle \frac{2}{N}}{\displaystyle \underset{j=0}{\overset{N}{}}}{}_{}{}^{\prime \prime }T_{j}^{}(\stackrel{~}{x}_k)T_j^{}(x).`$ (12)
and

$`{\displaystyle _1^x}f(t)dt`$ $``$ $`{\displaystyle \underset{k=0}{\overset{N}{}}}{}_{}{}^{\prime \prime }f(\stackrel{~}{x}_k){\displaystyle \frac{2}{N}}{\displaystyle \underset{j=0}{\overset{N}{}}}{}_{}{}^{\prime \prime }T_{j}^{}(\stackrel{~}{x}_k){\displaystyle _1^x}T_j(t)dt.`$ (13)
In matrix format, we have

$`\left[{\displaystyle _1^x}f(t)dt\right]`$ $``$ $`\stackrel{~}{S}\left[f\right],`$ (14)
$`\left[f^{}(x)\right]`$ $``$ $`\stackrel{~}{D}\left[f\right],`$ (15)
The elements of the column matrix $`\left[f\right]`$ are given by $`f(\stackrel{~}{x}_k),k=0(1)N`$. The right-hand side of Eqs. (14) and (15) give the values of the integral $`_1^xf(t)dt`$ and the derivative $`f^{}(x)`$ at the corresponding grid points, respectively. The actual values of the elements of the matrices $`\stackrel{~}{S}`$ and $`\stackrel{~}{D}`$ can be derived using Eqs. (12, 13).

In order to illustrate the Chebyshev algorithm, we consider again the case of a first-order nonlinear Volterra integro-differential equation of the form

$`𝐲^{}(x)=F[x,𝐲,Z[x;𝐲]],x[a,b],`$
$`Z[x;𝐲]={\displaystyle _a^x}K[x,t;𝐲(t)]dt,`$
with the initial condition

$`𝐲(a)=y_0.`$
Here we make no explicit restrictions on the actual form of the function $`F[x,𝐲,Z[x;𝐲]]`$, so both linear and nonlinear equations are included. We determine the unknown function $`𝐲(x)`$ using a perturbation approach: We start with an initial guess of the solution $`𝐲_0(x)`$ that satisfies the initial condition $`𝐲_0(a)=y_0`$, and write

$$𝐲(x)=𝐲_0(x)+ϵ(x),$$

with $`ϵ(x)`$ being a variation obeying the initial condition

$$ϵ(a)=0.$$
(16)
Hence, the original problem reduces to finding the perturbation $`ϵ(x)`$, and improving the initial guess in a iterative fashion.

We use the Taylor expansion of $`F[x,𝐲,Z[x;𝐲]]`$ about $`𝐲(x)=𝐲_0(x)`$ and keep only the linear terms in $`ϵ(x)`$ to obtain an equation for the variation $`ϵ(x)`$

$`ϵ^{}(x)`$ $``$ $`{\displaystyle \frac{F[x,𝐲,Z[x;𝐲]]}{𝐲(x)}}|_{𝐲(x)=𝐲_0(x)}ϵ(x)`$ (17)
$``$ $`{\displaystyle \frac{F[x,𝐲,Z[x;𝐲]]}{Z[x;𝐲]}}|_{𝐲(x)=𝐲_0(x)}{\displaystyle _a^x}{\displaystyle \frac{K[x,t;𝐲(t)]}{𝐲(x)}}|_{𝐲(x)=𝐲_0(x)}ϵ(t)\mathrm{d}t`$
$`=𝐲_0^{}(x)+F[x,𝐲_0(x),Z[x;𝐲_0(x)]].`$
Equation (17) is of the general form (18)

$$ϵ^{}(x)=q[x,ϵ(x)]+r(x),$$
(18)
where

$`q[x,ϵ(x)]={\displaystyle \frac{F[x,𝐲,Z[x;𝐲]]}{𝐲(x)}}|_{𝐲(x)=𝐲_0(x)}ϵ(x)`$
$`+{\displaystyle \frac{F[x,𝐲,Z[x;𝐲]]}{Z[x;𝐲]}}|_{𝐲(x)=𝐲_0(x)}{\displaystyle _a^x}{\displaystyle \frac{K[x,t;𝐲(t)]}{𝐲(x)}}|_{𝐲(x)=𝐲_0(x)}ϵ(t)\mathrm{d}t,`$
and

$`r(x)`$ $`=`$ $`𝐲_0^{}(x)F[x,𝐲_0(x),Z[x;𝐲]],`$
together with the initial condition given by (16). We replace Eqs. (18) and (16) by an integral equation, obtained by integrating Eq. (18) and using the initial condition (16) to choose the lower bound of the integral. We obtain

$`ϵ(x)={\displaystyle _a^x}q[t,ϵ(t)]dt+{\displaystyle _a^x}r(t)dt,`$ (19)
which is in fact a *linear* Volterra integral equation of the second kind. Using the techniques developed in the previous section to calculate integrals, the integral equation (19) can be transformed into a linear system of equations. A practical implementation of this algorithm is illustrated via a test problem in the following section.

## 4 Test problem

Following Shaw , we consider the test problem

$`𝐲(x)=xe^{1𝐲(x)}{\displaystyle \frac{1}{(1+x)^2}}x{\displaystyle _0^x}{\displaystyle \frac{x}{(1+t)^2}}e^{1𝐲(t)}dt,`$ (20)
$`𝐲(0)=y_0=1,x[0,1],`$ (21)
which has the exact solution

$`𝐲(x)={\displaystyle \frac{1}{1+x}}.`$ (22)
We shall use the initial guess $`𝐲_0(x)=y_0\mathrm{cos}(x)`$, so that $`𝐲_0(0)=y_0`$. The equation for the variation $`ϵ(x)`$ is

$`ϵ(x){\displaystyle _0^x}te^{1𝐲_0(t)}ϵ(t)dt+{\displaystyle _0^x}ds{\displaystyle _0^s}{\displaystyle \frac{se^{1𝐲_0(t)}}{(1+t)^2}}ϵ(t)dt=`$ (23)
$`𝐲_0(x)+y_0+{\displaystyle _0^x}\left[te^{1𝐲_0(t)}{\displaystyle \frac{1}{(1+t)^2}}t\right]dt{\displaystyle _0^x}ds{\displaystyle _0^s}{\displaystyle \frac{se^{1𝐲_0(t)}}{(1+t)^2}}dt.`$
In matrix format and using the Chebyshev expansion presented above, the variation $`ϵ(x)`$ will be obtained as the solution of linear system of equations

$$A\left[ϵ\right]=C,$$
(24)
with matrices $`A`$ and $`C`$ given as

$`A_{ij}`$ $`=`$ $`\delta _{ij}\stackrel{~}{S}_{ij}\left[te^{1𝐲_0(t)}\right]_j`$
$`+\stackrel{~}{S}_{ik}\stackrel{~}{x}_k\stackrel{~}{S}_{kj}\left[{\displaystyle \frac{e^{1𝐲_0(t)}}{(1+t)^2}}\right]_j,i,j=0(1)N,`$
$`C_i`$ $`=`$ $`[𝐲_0(t)]_i+y_0+\stackrel{~}{S}_{ik}\left[te^{1𝐲_0(t)}{\displaystyle \frac{1}{(1+t)^2}}t\right]_k`$
$`\stackrel{~}{S}_{ik}\stackrel{~}{x}_k\stackrel{~}{S}_k\mathrm{}\left[{\displaystyle \frac{e^{1𝐲_0(t)}}{(1+t)^2}}\right]_{\mathrm{}}.`$
From a computational point of view the computer time is spent initializing the matrix elements $`A_{ij}`$ and $`C_j`$ on one hand, and finding the solution of (24) on the other. On the first matter, the calculation decouples nicely, and once we have the vector $`[y_0]`$, we can calculate $`\{C_i,A_{ij},j=0(1)N\}`$ in parallel for $`i=0(1)N`$. The algorithm is as follows:

1. calculate $`[y_0]=[y_0]+[ϵ]`$ ;
2. broadcast $`[y_0]`$ ;
3. do $`i=0,N`$ :

1. master to slave: send $`i`$ ;
    2. slave: compute $`\{C_i,A_{ij},j=0(1)N\}`$ ;
    3. slave to master: return $`\{C_i,A_{ij},j=0(1)N\}`$.

Regarding the second step, i.e. solving the linear system of equations, the best choice is to use the machine specific subroutines, which generally outperform hand-coded solutions. When such subroutines are not available, as in the case of a Linux based PC cluster for instance, one can use one of the MPI implementations available on the market. We shall see that the efficiency of the equation solver is critical to the success of the parallel implementation of the Chebyshev-expansion approach. In order to illustrate this aspect we perform two calculations, first using a LU factorization algorithm, and secondly using an iterative biconjugate gradient algorithm. These are standard algorithms  for solving systems of linear equations, but their impact on the general efficiency of the approach is quite different.

### 4.1 Serial case

Figure 1 depicts the average CPU time required to complete the calculation for the various methods. Figure 2 illustrates the convergence of the two numerical methods. The spectral character of the method based on Chebyshev polynomials allows for an excellent representation of the solution for $`N>12`$. We base our findings on a $`\sigma <10^{10}`$ criteria, where $`\sigma `$ denotes the sum of all absolute departures of the calculated values from the exact ones, at the grid points.

The number of iterations required to achieve the desired accuracy in the Chebyshev case is depicted in Fig. 3. The number of iterations becomes flat for $`N>12`$, and stays constant (17 iterations) even for very large values of N. The higher number of iterations corresponding to the lower values of N, represents an indication of a insufficient number of Chebyshev grid points: the exact solution cannot be accurately represented as polynomial of degree N for $`x[0,1]`$. It is interesting to note that for $`N=1216`$, a reasonable lower domain for the representation of the solution using Chebyshev polynomials, the reported CPU time is so small that for our test problem there is no real justification for porting the algorithm to a MPP machine. This situation will change for multi-dimensional problems such as those encountered in our nonequilibrium quantum field theory studies.

### 4.2 Parallel case

The LU factorization algorithm is an algorithm of order $`N^3`$ and consequently, most of the CPU time is spent solving the linear system of equations (see Fig. 4). As a consequence, a parallel implementation of the LU algorithm is very difficult. Figure 5 shows how the average CPU time changes with the available number of processors. Here we use a very simple MPI implementation of the LU algorithm as presented in reference . Even though we could certainly achieve better performance by employing a sophisticated LU equation solver, the results are typical. Since the actual size of the matrices involved is small, the communication overhead is overwhelming and the execution time does not scale with the number of processors.

Fortunately, even for dense matrices and small values of the number of grid points $`N`$, one can achieve a good parallel efficiency. By employing an iterative method such as the iterative biconjugate gradient method, one can render the time required to solve the system of linear equations negligible compared with the time required to initialize the relevant matrices, which in turn is only slightly more expensive than the initialization process of the LU factorization algorithm. The initialization process can be parallelized using the algorithm presented above and the results are depicted in Fig. 6.

It appears that by using the biconjugate gradient method the efficiency of the parallel code has improved considerably. However, the average CPU time saturates to give an overall speedup of 3.5 . This can be understood by analyzing the computation and communication requirements for our particular problem. The calculation cost to initialize the matrices $`A`$ and $`C`$ is roughly given by the number of floating-point multiplications and additions $`(7N^2+3N)T_{\text{calc}}`$, while the communication cost is given by $`(N^2+2N)T_{\text{comm}}`$. Therefore, the ratio of communication to computation is

$$\frac{N^2+2N}{7N^2+3N}\frac{T_{\text{comm}}}{T_{\text{calc}}}.$$

As in the finite-difference case, this ratio approaches a *constant* value as $`N`$ gets larger and it becomes apparent that the communication overhead is still a problem.

However, multi-dimensional applications such as those presented in  require complicated matrix element calculation. In such cases, the process of initializing the matrices $`A`$ and $`C`$ is quite involved, and the ratio of the communication time relative to the computation time becomes favorable. In addition, the matrix $`A`$ becomes sparse and the size of the linear system of equations is substantially larger, thus one can also take advantage of existing parallel implementation of the iterative biconjugate gradient algorithm . Such problems benefit heavily from an adequate parallelization of the code. We will discuss such an example in the following section.

## 5 Volterra-like integral equations for a two-point Green function

Schwinger, Bakshi, Mahanthappa, and Keldysh  have established how to formulate an initial value problem in quantum field theory. The formalism is based on a generating functional, and the evolution of the density matrix requires both a forward evolution from zero to $`t`$ and a backward one from $`t`$ to zero. This involves  both positive and negative time ordered operators in the evolution of the observable operators and the introduction of two currents into the path integral for the generating functional. Time integrals are then replaced by integrals along the closed time path (CTP) in the complex time plane shown in Fig. 7. We have

$$_𝒞F(t)dt=_{0:𝒞_+}^{\mathrm{}}F_+(t)dt_{0:𝒞_{}}^{\mathrm{}}F_{}(t)dt.$$
(25)
Using the CTP contour, the full closed time path Green function for the two point functions is:

$$𝒢(t,t^{})=𝒢_>(t,t^{})\mathrm{\Theta }_𝒞(t,t^{})+𝒢_<(t,t^{})\mathrm{\Theta }_𝒞(t^{},t),$$

in terms of the Wightman functions, $`𝒢_{>,<}(t,t^{})`$, where the CTP step function $`\mathrm{\Theta }_𝒞(t,t^{})`$ is defined by:

$$\mathrm{\Theta }_𝒞(t,t^{})=\{\begin{array}{cc}\mathrm{\Theta }(t,t^{})\hfill & \text{for }t\text{ on }𝒞_+\text{ and }t^{}\text{ on }𝒞_+,\hfill \\ 0\hfill & \text{for }t\text{ on }𝒞_+\text{ and }t^{}\text{ on }𝒞_{},\hfill \\ 1\hfill & \text{for }t\text{ on }𝒞_{}\text{ and }t^{}\text{ on }𝒞_+,\hfill \\ \mathrm{\Theta }(t^{},t)\hfill & \text{for }t\text{ on }𝒞_{}\text{ and }t^{}\text{ on }𝒞_{}.\hfill \end{array}$$
(26)
For complete details of this formalism and various applications, we refer the reader to the original literature , and we confine ourselves to discussing how our Chebyshev-expansion approach is applied to the computation of the two-point Green function.

For simplicity we consider now the quantum mechanical limit of quantum field theory (0+1 dimensions). In this limit, we are generally faced with the problem of numerically finding the solution of equation

$$𝒢(t,t^{})=G(t,t^{})_𝒞dt^{\prime \prime }Q(t,t^{\prime \prime })𝒢(t^{\prime \prime },t^{}),$$
(27)
Here, the Green functions, $`𝒢(t,t^{})`$ and $`G(t,t^{})`$, are symmetric in the sense that $`𝒢_>(t,t^{})=𝒢_<(t^{},t)`$, and obey the additional condition

$$𝒢_{>,<}(t,t^{})=𝒢_{<,>}^{}(t,t^{})=𝒢_{<,>}(t^{},t).$$
(28)
The function $`Q(t,t^{})`$ obeys less stringent symmetries

$$Q_{>,<}(t,t^{})=Q_{<,>}^{}(t,t^{})Q_{<,>}(t^{},t),$$
(29)
which is always the case when $`Q(t,t^{})`$ has the form

$$Q(t,t^{})=_𝒞dt^{\prime \prime }A(t,t^{\prime \prime })B(t^{\prime \prime },t^{}),$$
(30)
where $`A(t,t^{})`$ and $`B(t,t^{})`$ satisfy (28).

We can further write Eq. (28) as

$`e\{𝒢_>(t,t^{})\}`$ $`=`$ $`e\{𝒢_<(t,t^{})\},`$ (31)
$`m\{𝒢_>(t,t^{})\}`$ $`=`$ $`m\{𝒢_<(t,t^{})\},`$ (32)
or

$`𝒢_>(t,t^{})𝒢_<^{}(t,t^{})`$ $`=`$ $`2e\{𝒢_>(t,t^{})\},`$ (33)
$`𝒢_>(t,t^{})+𝒢_<^{}(t,t^{})`$ $`=`$ $`2m\{𝒢_>(t,t^{})\}.`$ (34)
Hence, a Green function $`𝒢(t,t^{})`$ is fully determined by the component $`𝒢_>(t,t^{})=e\{𝒢_>(t,t^{})\}+im\{𝒢_>(t,t^{})\}`$, with $`t^{}t`$. Thus, in order to obtain the solution of Eq. (27), we only need to solve

$`𝒢_>(t,t^{})`$ $`=`$ $`G_>(t,t^{})2{\displaystyle _0^t}dt^{\prime \prime }e\{Q_>(t,t^{\prime \prime })\}𝒢_>(t^{\prime \prime },t^{})`$ (35)
$`+2{\displaystyle _0^t^{}}dt^{\prime \prime }Q_>(t,t^{\prime \prime })e\{𝒢_>(t^{\prime \prime },t^{})\}.`$
We separate the real and the imaginary part of (35) and obtain the system of integral equations

$`e\{𝒢_>(t,t^{})\}`$ $`=`$ $`e\{G_>(t,t^{})\}2{\displaystyle _0^t}dt^{\prime \prime }e\{Q_>(t,t^{\prime \prime })\}e\{𝒢_>(t^{\prime \prime },t^{})\}`$ (36)
$`+2{\displaystyle _0^t^{}}dt^{\prime \prime }e\{Q_>(t,t^{\prime \prime })\}e\{𝒢_>(t^{\prime \prime },t^{})\}`$
$`m\{𝒢_>(t,t^{})\}`$ $`=`$ $`m\{G_>(t,t^{})\}2{\displaystyle _0^t}dt^{\prime \prime }e\{Q_>(t,t^{\prime \prime })\}m\{𝒢_>(t^{\prime \prime },t^{})\}`$ (37)
$`+2{\displaystyle _0^t^{}}dt^{\prime \prime }m\{Q_>(t,t^{\prime \prime })\}e\{𝒢_>(t^{\prime \prime },t^{})\}.`$
The above system of equations must be solved for $`t^{}t`$. The two equations are independent, which allows us to solve first for the real part of $`𝒢_>(t,t^{})`$, and then use this result to derive the imaginary part of $`𝒢_>(t,t^{})`$.

Despite their somewhat unusual form, the above equations are two-dimensional Volterra-like integral equations and our general discussion regarding the Chebyshev spectral method applies. We will perform a multi-step implementation of the formalism. Let

$$t_i=t_{i_0(N1)+i_1},1i_1N,$$

be the grid location corresponding to the collocation point $`i_1`$ of the interval labelled $`i_0+1`$. Then, the discrete correspondent of Eq. (35) is

$`𝒢_>(t_i,t_j)`$ $`=`$ $`G_>(t_i,t_j)`$
$`{\displaystyle \underset{k_0=0}{\overset{i_01}{}}}{\displaystyle \underset{k_1=1}{\overset{N}{}}}[2\stackrel{~}{S}_{Nk_1}]e\{Q_>(t_i,t_{k[=k_0(N1)+k_1]})\}𝒢_>(t_k,t_j)`$
$`{\displaystyle \underset{k_1=1}{\overset{N}{}}}[2\stackrel{~}{S}_{i_1k_1}]e\{Q_>(t_i,t_{k[=i_0(N1)+k_1]})\}𝒢_>(t_k,t_j)`$
$`+{\displaystyle \underset{k_0=0}{\overset{j_01}{}}}{\displaystyle \underset{k_1=1}{\overset{N}{}}}[2\stackrel{~}{S}_{Nk_1}]Q_>(t_i,t_{k[=k_0(N1)+k_1]})e\{𝒢_>(t_k,t_j)\}`$
$`+{\displaystyle \underset{k_1=1}{\overset{N}{}}}[2\stackrel{~}{S}_{j_1k_1}]Q_>(t_i,t_{k[=j_0(N1)+k_1]})e\{𝒢_>(t_k,t_j)\},`$
with $`t_jt_i`$.

We will refer now to Figs. 8 and 9. Equation (5) involves values of $`𝒢_>(t_k,t_j)`$, for which $`t_j>t_k`$. In such cases, we use the symmetry $`𝒢_>^{}(t_j,t_k)`$, which relates to the values the two-point function located in the domain of interest. For the time interval $`(i_0+1)`$ the size of the linear system of equations we need to solve is

$`{\displaystyle \frac{1}{2}}(i_0+1)(N1)[(i_0+1)(N1)+1]{\displaystyle \frac{1}{2}}i_0(N1)[i_0(N1)+1]`$
$`=i_0(N1)^2+{\displaystyle \frac{1}{2}}N(N1),`$
or of order $`(i_0+\frac{1}{2})(N1)^2`$. In practice, the value of $`N`$ is taken between 16 and 32.

Tables 1 and 2 summarize the number of floating-point operations performed in order to compute the non-vanishing matrix elements corresponding to a given $`i`$ and $`j,(j<i)`$.

We can now calculate the ratio of communication to computation time, by noticing that the numbers in the tables above get multiplied by N, corresponding to the number of collocation points in each time step and summing over the number of steps, i.e. we evaluate

$$N\left[\text{if}j>i_0(N1)\right]+N\underset{j_0=1}{\overset{i_0}{}}\left[\text{if}ji_0(N1)\right].$$

In Table 3 we summarize all relevant estimates regarding the computation cost for a fixed value of $`i`$. In order to estimate the *total* communication and computation cost, respectively, these numbers must be multiplied by an additional factor of $`N`$, corresponding to the number of possible values of $`i`$ in a time step. This factor is not relevant for estimating the communication overhead, but it must be remembered when one infers the sparsity of the corresponding system of equations.

To conclude we observe that the communication to computation ratio approaches

$$\frac{1}{2(i_0+1)}\frac{T_{\text{comm}}}{T_{\text{calc}}}$$

for large values of $`i_0`$. Therefore for this problem the communication overhead is reduced substantially in the later stages of the calculation. In practice, this ratio is actually much better, as we compute the functions $`G(t,t^{})`$ and $`Q(t,t^{})`$ on the fly, and this adds considerably to the computational effort. Finally the sparsity of the resulting systems of equations goes to $`2/(i_0N)`$ for large values of $`i_0`$ and $`N`$, which supports our choice for an iterative equation solver.

## 6 Conclusions

We have presented a numerical method suitable for solving non-linear integral and integro-differential equations on a massively multiprocessor machine. Our approach is essentially a standard perturbative approach, where one calculates corrections to an initial guess of the solution. The initial guess is designed to satisfy the boundary conditions, and corrections are expanded out in a complete basis of N Chebyshev polynomials on the grid of (N+1) extrema of $`T_N(x)`$, the Chebyshev polynomial of first kind of degree N. The spectral character of the convergence of the Chebyshev-expansion approach is the key element in keeping low the number of grid points. From a computational point of view, each iteration involves two stages, namely initializing the relevant matrices and solving the linear system of equations. Both stages can be rendered parallel in a suitable manner, and the efficiency of the code increases when applied to complicated multi-step, multi-dimensional problems.

The algorithm discussed in this paper represents the backbone of current investigations of the equilibrium and nonequilibrium properties of various phenomenological Lagrangeians. In particular we are interested in studying the properties of the chiral phase transition at finite density for a 2+1 dimensional four-fermion interaction as well as the dynamics od 2-dimensional QCD, with the ultimate goal of indirectly obtaining insights regarding the time evolution of a quark-gluon plasma produced following a relativistic heavy-ion collision.

The work of B.M. was supported in part by the U.S. Department of Energy, Nuclear Physics Division, under contract No. W-31-109-ENG-38. The work of R.S. was supported in part by the Natural Sciences and Engineering Research Council of Canada under grant No. OGP0170170. Parallel calculations are made possible by grants of time on the parallel computers of the Mathematics and Computer Science Division, Argonne National Laboratory. B.M. would like to acknowledge useful discussions with John Dawson and Fred Cooper.

## References


format: text

text:
# UHECR as Decay Products of Heavy Relics? The Lifetime Problem

## 1 Introduction

Ultrahigh Energy Cosmic Rays (UHECR) have been observed by a number of experiments at energies above $`10^{20}`$ eV. The forthcoming cosmic rays detectors as the Auger array, the EUSO and OWL space observatories are expected to greatly improve our present knowledge on the UHECR gathered from Fly’s Eye, HiRes, AGASA and previous detectors.

Top-down scenarii for UHECR are based on heavy relics from the early universe which are assumed to decay at the present time. In all cases, whatever the nature of the objects: heavy particles, topological and non-topological solitons, black-holes, microscopic fundamental strings, cosmic defects etc., one has to fine tune the lifetime of these objects to be the age of the universe.

We provide an unified description for the quantum decay formula of unstable particles which encompass all the above mentioned cases, as well as the particle decays in the standard model (muons, Higgs, etc). In all cases the decay rate can be written as,

$$\mathrm{\Gamma }=\frac{g^2m}{\text{numerical factor}}$$
(1)
where $`g`$ is the coupling constant, $`m`$ is the typical mass in the theory (it could be the mass of the unstable particle) and the numerical factor contains often relevant mass ratios for the decay process.

The key drawback of all top-down scenarii is the lifetime problem. The ad-hoc requirement of a lifetime of the order the age of the universe for the heavy particles implies an operator with a very high dimension describing the decay, and/or an extremely small coupling constant.

Heavy relics could have been formed by the end of inflation at typical GUT’s energy scales, but their natural lifetime would be of the order of microscopic times typically associated to GUT’s energy scales .

UHECR may result from the acceleration of protons and ions by shock-waves in astrophysical plasmas (Fermi acceleration mechanism). Large enough sources can accelerate particles to the energies of the observed UHECR. Sources in the vicinity of our galaxy as hot spots of radio galaxies (working surfaces of jets and the inter galactic medium) and blazars (active galactic nuclei with relativistic jet directed along the line of sight) as BL Lacertae can evade the GZK bound.

## 2 Topological solitons, non-topological solitons and heavy particles

Stable solutions in classical field theory (as monopoles) become (heavy) particles in quantum field theory. There is no difference at the quantum level between heavy particles associated to a local field and those associated to classical stable solutions.

The stability of classical solutions in field theory is a highly nontrivial issue. There are basically two types of solutions: topological and non-topological. Topological classical solutions have associated a non-zero topological number (topological charge) which vanishes for the vacuum. If there is a lower bound for the energy of the solution involving this topological number the classical solution is stable. This is the case for kinks in one space dimension scalar theories, vortices in the two-dimensional Higgs model, monopoles in the three dimensional Georgi-Glashow model, Hopf solitons in appropriate three dimensional scalar models. In all known cases, classical stability comes together with quantum stability.

Gravitational analogues of these classical solutions exist in the euclidean (imaginary time) regime : they are black-holes in three space dimensions (with periodicity in the imaginary time), which are gravitational analogues of electric type monopoles and Taub-Nut’s in four space dimensions (gravitational analogues of magnetic type monopoles). The topological charges here are related to the temperature and magnetic charge of the solutions, respectively.

It must be stressed that the mere presence of a conserved topological number does not guarantee the stability of the corresponding classical solution. The energy must be related with the topological number in question such that a non-zero topological number implies a non-zero energy. Otherwise, a classical solution possessing non-zero topological number can decay into lighter particles.

In other words, the topological charge may be disconnected from the dynamics and it can decay in the course of the evolution. A topological soliton may collapse loosing its topological charge. This does not happen when the topological charge bounds the energy from below.

Non-topological solitons are stable thanks to a conserved $`U(1)`$ charge of ‘electric’ type. Again, the mere presence of a conserved $`U(1)`$ charge does not guarantee stability for charged particles except for the lightest one. Let us call $`m`$ the mass of the lightest charged particle and let us take its $`U(1)`$ charge as unit of charge. Assume that there are heavier particles with mass $`M>m`$ and charge $`Q>1`$ with $`M=M(Q)`$. A sufficient condition for quantum stability is

$$M(Q)<mQ,$$

since a particle with charge $`Q`$ and mass larger or equal than $`mQ`$ can always decay into $`Q`$ particles of mass $`m`$ and unit charge respecting charge and energy conservation.

It must be stressed that in quantum theory all non-forbidden process do happen.

## 3 Quantum Decay of Heavy Particles

Typically, the decay of a heavy particle with mass $`m_X`$ can be described by an effective interaction lagrangian formed by the local field $`X(x)`$ associated to this heavy particle times the lighter fields in which it decays. Let us take the muon decay which is a well known case. Notice the mass of the muon $`m_\mu =206.8m_em_e`$.

The effective Fermi lagrangian can be written as

$$_I=\frac{G_F}{\sqrt{2}}\overline{\psi }_{\nu _\mu }\gamma ^\alpha (1+\gamma _5)\psi _\mu \overline{\psi }_{\nu _e}\gamma _\alpha (1+\gamma _5)\psi _{\nu _e}$$
(2)
where $`\psi _\mu `$ stands for the muon field, $`\psi _{\nu _e}`$ and $`\psi _{\nu _\mu }`$ for the electron neutrino and muon neutrino fields, respectively. The Fermi coupling $`G_F`$ has the dimension of an inverse square mass.

The muon width $`\mathrm{\Gamma }_\mu `$ describing the decay is then given by

$$\mathrm{\Gamma }_\mu =\frac{G_F^2m_\mu ^5}{192\pi ^3}$$

The Fermi coupling can be related to the W-mass as follows

$$\frac{G_F}{\sqrt{2}}=\frac{g^2}{8m_W^2}$$

where $`g`$, the standard model coupling, is dimensionless. Thus,

$$\mathrm{\Gamma }_\mu =\frac{g^4m_\mu }{6144\pi ^3}\left(\frac{m_\mu }{m_W}\right)^4$$
(3)
As we shall see below, Eq.(3) has the generic structure for the decay width of an unstable particle.

For the muon decay, the monomial interaction in the effective lagrangian (2) has dimension six in mass units.

An analogous example is the Higgs decay into muons, neutrinos, $`W^\pm `$ and the $`Z^0`$. Notice that the Higgs mass $`m_H`$ must be higher than the $`W^\pm `$ mass $`m_W`$ and the $`Z^0`$ mass. The lagrangian as given by the standard model is here

$$2g\mathrm{sin}\theta _WM_WHW_\mu ^+W_{}^\mu $$
(4)
and a similar expression for the coupling with the $`Z`$. Here $`\theta _W`$ stands for Weinberg’s angle.

One finds for the Higgs decay rate,

$$\mathrm{\Gamma }_{Higgs}=\frac{3g^2}{128\pi }m_H\left(\frac{m_H}{m_W}\right)^2$$
(5)
where we consider for simplicity the case $`M_HM_W`$. In this case the monomial interaction in the effective lagrangian (4) has dimension three in mass units.

Notice that in both cases, eq.(3) and eq.(5), the width grows as a positive power of the mass of the decaying particle.

Let us consider an effective lagrangian containing a local monomial of dimension $`n`$ (in mass units)

$$_I=\frac{g}{M^{n4}}X\mathrm{\Theta }.$$
(6)
Here the field $`X`$ is associated to the decaying particle of mass $`m_X`$ and $`\mathrm{\Theta }`$ stands for the product of fields coupled to it.

Then, the decay rate for a particle of mass $`m_X`$ takes the form

$$\mathrm{\Gamma }=\frac{g^2}{\text{numerical factor}}m_X\left(\frac{m_X}{M}\right)^{|2n8|}$$
(7)
$`\mathrm{\Gamma }_\mu `$ eq.(3) and $`\mathrm{\Gamma }_{Higgs}`$ eq.(5) correspond to $`n=6`$ and $`n=3`$, respectively.

## 4 Quantum Decay of Solitons

The mass of classical soliton solutions (as magnetic monopoles in unified theories) are of the form

$$M_{sol}=\frac{\mu }{g^2}$$

where $`\mu `$ is the mass of the basic fields in the lagrangian and $`g`$ their dimensionless coupling. For small coupling these objects are much heavier than the particles associated to the basic fields in the lagrangian.

Quantum mechanically the soliton mass acquires corrections of order $`g^0`$ and higher. To one-loop level one finds

$$M_{sol}=\frac{\mu }{g^2}+\frac{1}{2}\underset{n}{}\left[\omega _n\omega _n^0\right],$$
(8)
where $`\omega _n`$ stands for the frequency of oscillations around the soliton. These oscillations are close but not identical to the frequency of oscillations around the vacuum $`\omega _n^0`$. The sum in eq.(8) yields a finite result proportional to $`\mu `$.

Now, if the classical solution is unstable some of the frequencies $`\omega _n`$ develops an imaginary part $`i\mu \beta `$ where $`\beta `$ is a pure number. Hence,

$$ImM_{sol}=\beta \mu \text{and}ReM_{sol}=\frac{\mu }{g^2}+𝒪(g^0)$$
(9)
and

$$\mathrm{\Gamma }_{sol}=ImM_{sol}=g^2\beta ReM_{sol}$$
(10)
We see that the width $`\mathrm{\Gamma }_{sol}`$ has a similar structure than for heavy particles in the previous section.

The term $`𝒪(g^0)`$ in eq.(9) stand for the first quantum correction to the mass. Notice that we choose $`\mathrm{}=1`$ which is absorbed in $`g^2`$.

## 5 Quantum Decay of Fundamental Strings

The decay of closed strings in string theory has been computed to the dominant order (one string loop). Assuming the closed string in an $`N`$th excited state, it can decay into lower excited states including the graviton and the dilaton. The mass of this quantum string is given by

$$m^2=32\pi TN$$

where $`T`$ is the string tension $`T=1/(4\pi \alpha ^{})`$ and $`\alpha ^{}`$ the string constant. The length of such string is given by $`L=2\alpha ^{}m`$.

One finds for the total width for string decay,

$$\mathrm{\Gamma }_{string}=\frac{\kappa ^2\sqrt{TN}}{\text{numerical factor}}$$
(11)
where the dimensionless coupling $`\kappa `$ is given by,

$$\kappa =48\pi \sqrt{2GT}$$

The total width can be then rewritten as,

$$\mathrm{\Gamma }_{string}=\frac{\kappa ^2m}{1083\times \text{numerical factor}}.$$
(12)
This formula again has the same structure as the previous widths (7) and (10) once we identify $`g=\kappa ,m_X=ReM_S=m`$.

Eq.(12) can be rewritten as,

$$\mathrm{\Gamma }_{string}=42\frac{GTm}{\text{numerical factor}}=\frac{21}{16\pi }\frac{Gm^3}{\text{numerical factor}}=\frac{21\sqrt{2}}{\text{numerical factor}}\frac{\sqrt{N}G}{\alpha ^{3/2}}$$
(13)
## 6 Quantum Decay of Black Holes

As it is known, in the context of field theory black holes decay semiclassically through thermal emission at the Hawking temperature

$$T_{BH}=\frac{\mathrm{}c}{4\pi k_B}\frac{1}{R_s},R_s=\frac{2GM}{c^2}$$

($`M`$ being the black hole mass and $`k_B`$ the Boltzmann constant).

Black hole emission follows a ‘gray body’ spectrum (the ‘filter’ being the black hole absorption cross section $`R_s^2`$). The mass loss rate in this process can be estimated following a Stefan-Boltzmann relation,

$$\frac{dM}{dt}=\sigma R_s^2T_{BH}^4T_{BH}^2$$

where $`\sigma `$ is a constant. Thus, the black hole decay rate is

$$\mathrm{\Gamma }_{BH}=\left|\frac{1}{M}\frac{dM}{dt}\right|GT_{BH}^3\frac{G}{R_s^3}$$

As evaporation proceeds, the black hole temperature increases till it reaches the string temperature

$$T_{string}=\frac{\mathrm{}c}{k_B}\frac{1}{bL_s},L_s=\sqrt{\frac{\mathrm{}\alpha ^{}}{c}}$$

($`L_s`$ being the fundamental string length and $`b`$ a constant exclusively depending on the spacetime dimensionality and the string model chosen.) The black hole enters its string regime $`T_{BH}T_{string},R_sL_s`$, becomes a string state and decays with a width

$$\mathrm{\Gamma }_{BH}GT_s^3\frac{G}{\alpha ^{3/2}}\mathrm{\Gamma }_{string}.$$

Notice that this formula is similar to eq.(13) and again has the generic structure of the widths eq.(7)-(10) and (12) if one identifies $`g=\kappa ,m_X=ReM_S=m`$.

We consider here both fundamental strings and black holes since their decay rates can be nicely recasted as in eq.(1) independently of whether or not they may be considered as candidate sources of UHECR.

## 7 Particles Lifetime and the Age of the Universe

Heavy particles with masses in the GUT scale can be produced in large numbers during inflation and just after inflation. The production mechanism is parametric or spinodal amplification in the inflaton field. That is, linear resonance of the quantum modes of the heavy field in the background or condensate of the inflaton. In addition, non-linear quantum phenomena play a crucial role and can enhance the particle production. Such non-linear production is of the same order of magnitude as the gravitational production of particles by the time dependent metric.

Once these heavy particles are produced they must have a lifetime of the order of the age of the universe in order to survive in the present universe and decay into UHE cosmic rays. Only in the early universe the production of such heavy objects is feasible due to their large mass.

Moreover, in order to be the source of UHECR, these particles must have a mass of the observed UHECR, namely $`m_X>10^{21}`$ eV = $`10^{12}`$ GeV.

Let us assume that the effective lagrangian (6) describes the decay of the $`X`$ particles. Their lifetime will be given by eq.(7)

$$\tau _X=\frac{\text{numerical factor}}{g^2}\frac{1}{m_X}\left(\frac{M}{m_X}\right)^{2n8}=\frac{\text{numerical factor}}{g^2}\frac{1}{m_X}\mathrm{\hspace{0.33em}10}^{6(n4)}$$

where we set a GUT mass $`M=10^{15}`$ GeV. The age of the universe is $`\tau _{universe}\mathrm{2\hspace{0.33em}10}^{10}`$years and we have to require that $`\tau _X>\tau _{universe}`$. Therefore,

$$10^{54}<\frac{\text{numerical factor}}{g^2}\mathrm{\hspace{0.33em}10}^{6(n4)}\text{or}\mathrm{log}_{10}g<3(n13)$$
(14)
and we dropped the numerical factor in the last step.

For $`g1`$, eq.(14) requires an operator $`\mathrm{\Theta }`$ with dimension at least thirteen in the effective lagrangian (6) which is a pretty high dimension.

That is, one needs to exclude all operators of dimension lower than thirteen in order to extremely suppress the decay. Clearly, one may accept lower dimension operators $`\mathrm{\Theta }`$ paying the price of a small coupling $`g`$. For example: $`g=10^9`$ and $`n=10`$ fullfil the above bound still being a pretty high dimension operator. Notice that a moderate $`n`$ as $`n=4`$ lowers the coupling to $`g10^{27}`$.

In summary, a heavy X-particle can survive from the early universe till the present times if one chooses

* an extremely small coupling $`g`$

and/or
* an operator $`\mathrm{\Theta }`$ with high enough dimension

None of these assumptions can be supported by arguments other than imposing a lifetime of the age of the universe to the X-particle. That is, the lifetime must be here fine tunned. That is, one has to built an ad-hoc lagrangian to describe the X-particle decay. Indeed, a variety of ad-hoc lagrangians have been proposed in the literature together with the symmetries which can adjust a wide variety of lifetimes.

It must be recalled that no known (weakly broken) symmetry protects the X-particle from decaying rapidly, except for supersymmetry. However, if supersymmetry would be invoked in this context, that would imply that the supersymmetry scale is at the GUT scale or beyond. It must be also noticed that the natural lifetime for particles of such a mass is the GUT scale, that is typically $`10^{28}`$ sec. - $`10^{35}`$ sec.

## 8 Cosmic Defects and Heavy Particles

Closed vortices from abelian and non-abelian gauge theories are not topologically stable in $`3+1`$ space-time dimensions. Static vortices in $`3+1`$ space-time dimensions just collapse to a point since their energy is proportional to their length. They do that in a very short (microscopic) time.

It must be noticed that only a restricted set of spontaneously broken non-abelian gauge theories exhibit vortex solutions. For example, there are no topologically stable vortices in the standard $`G=SU(3)\times SU(2)\times U(1)`$ model in $`3+1`$ space-time dimensions just because $`\mathrm{\Pi }_1(G)`$ and $`\mathrm{\Pi }_2(G)`$ are trivial for such group manifolds. \[For a recent review see \]. Grand unified theories may or may not posses vortex solutions in $`2+1`$ space-time dimensions depending under which representations of the gauge group belong the Higgs fields.

The existence of cosmic string networks is not established although they have been the subject of many works. In case such networks would have existed in the early universe they may have produced heavy particles X of the type discussed before and all the discussion on their lifetime applies here. The discussion on the lifetime problem also applies to rotating superconducting strings which have been proposed as classically stable objects.

Cosmic strings are closed vortices of horizon size. In $`3+1`$ space-time dimensions, strings collapse very fast except if they have horizon size in which case their lifetime would be of the order of the age of the universe. However, such horizon size cosmic strings are excluded by the CMB anisotropy observations and by the isotropy of cosmic rays. Such gigantic objects behaves classically whereas microscopic closed strings (for energies $`<M_{Planck}=10^{19}`$ GeV) behave quantum mechanically.

In summary, a key point here is the unstability of topological defects in $`3+1`$ space-time dimensions. Unless one chooses very specific models , topological defects decay even classically with a short lifetime. They collapse to a point at a speed of the order of the speed of light in $`3+1`$ space-time dimensions.

## 9 Decay Products of Heavy Particles and Final Conclusions

In summary, if the X-particles, whatever their origin and type could be made sufficiently stable to survive till now, then their decay products could provide the UHECR observed today. However, the X-particles lifetime of the order of the age of the universe must be imposed ad-hoc i. e. fine tuned while the natural lifetime for those particles should be extremely short about $`10^{28}`$ sec at most.

Various GUTs contain candidates for X-particles of masses around the GUT scale ranging approximately from $`10^{12}`$ Gev to $`10^{16}`$ Gev depending on the model. These particles could have been produced naturally in the early universe typically by the end of inflation. Analogously, topological defects, fundamental strings and primordial black holes could have been formed in the early universe. The hard job, however, is to have these heavy objects still present and decaying today. Instead of that, it seems more natural that the X-particles and the other heavy objects above mentioned decayed in the early universe shortly after being formed, having lifetimes corresponding to their respective energy scales. Their decay products will then form relic primordial backgrounds as graviton, neutrino and dilaton backgrounds, as we know now the relic photon CMB background. Those backgrounds could have characteristic detectable spectra and signatures containing informations about the early universe.


text:
# Discovery of complex narrow X–ray absorption features from the low-mass X–ray binary GX 13+1 with XMM-Newton

## 1 Introduction

The bright and persistent low-mass X–ray binary (LMXRB) GX 13$`+`$1 (4U 1811–17) is an X–ray burst source (Fleischman f:85 (1985); Matsuba m:95 (1995)) with radio (Grindlay & Seaquist gs:86 (1986)) and infrared counterparts (Naylor et al. ncl:91 (1991); Charles & Naylor cn:92 (1992)). It displays a 15% modulation in X–ray intensity with a period of 24.7$`\pm `$1 days. The X–ray intensity is anti-correlated with the hardness ratio (Corbet c:96 (1996)). Recent estimates of the spectral type of the companion (an evolved late-type K5 iii star; Bandyopadhyay et al. b:99 (1999)), together with the binary system parameters are consistent with this modulation being the orbital period, although further confirmation is required (Bandyopadhyay et al. b:02 (2002)). GX 13$`+`$1 has been classified as an atoll source, but in some respects its properties are more similar to those of the Z-sources (Homan et al. h:98 (1998); Schulz et al. s:89 (1989); Bandyopadhyay et al. b:99 (1999)), having persistent radio emission at a similar level to the Z-sources (Fender & Hendry fh:00 (2000)), and a high accretion rate of $``$0.1 Eddington. A quasi-periodic oscillation (QPO) was discovered at 57–69 Hz by Homan et al. (h:98 (1998)). Distance estimates for GX 13$`+`$1 range from 4.7 to 7.0 kpc (Christian & Swank cs:97 (1997)). Infrared spectroscopy favors a distance of 7$`\pm `$1 kpc (Bandyopadhyay et al. b:99 (1999)).

A spectral survey of Fe-K emission from 20 LMXRB with ASCA is reported in Asai et al. (a:00 (2000)). Significant iron emission was detected from GX 13$`+`$1 with a line centroid of $`6.44\pm 0.05`$ keV, an equivalent width (EW) of $`37_8^{+10}`$ eV and a line width $`\sigma `$$``$130$`\pm `$40 eV. The X–ray continuum could be modeled with a two-component model, consisting of a multicolor disk-blackbody with a temperature, $`kT`$, at the inner disk radius of 0.76 keV and a blackbody with a $`kT`$ of 1.34 keV, absorbed by a column density, $`N_\mathrm{H}`$, of $`\mathrm{2.9\hspace{0.17em}10}^{22}`$ atom cm-2. The source flux (1–10 keV, corrected for absorption) was 1.2 10-8 erg cm-2 s-1, which translates into a luminosity of $`7`$1037 erg s-1 (at 7 kpc). Recently, Ueda et al. (u:01 (2001)) re-analyzed the ASCA Solid-State Imaging Spectrometer (SIS) observation of GX 13$`+`$1 and discovered that the spectral region around Fe is more complex than had been realized. Their fits required the presence of a narrow ($`\sigma <70`$ eV) absorption line at $`7.01\pm 0.03`$ keV with an EW of 35$`\pm `$8 eV, which they interpreted as resonant scattering of the K$`\alpha `$ line from Fe xxvi ions with a column density $`\text{ }>`$$`10^{18}`$ cm-2, a narrow ($`\sigma <220`$ eV) emission line at $`6.42\pm 0.08`$ keV with an EW of $`19\pm 8`$ eV and an absorption edge at $`7.61\pm 0.13`$ keV with an optical depth of $`0.13\pm 0.05`$. This was the first detection of an Fe absorption line in an X–ray binary known to contain a neutron star. Previously, the only X–ray binaries known to exhibit such Fe absorption features were two superluminal jet sources GRO J1655-40 (Ueda et al. u:98 (1998); Yamaoka et al. y:01 (2001)) and GRS 1915+105 (Kotani et al. k:00 (2000); Lee et al. l:01 (2001)). As Ueda et al. (u:01 (2001)) comment, the interpretation of the edge is not simple and it could be a combination of an absorption edge from moderately ionized Fe together with line absorption from a number of prominent K$`\alpha `$ or K$`\beta `$ lines of Fe and Ni.

XMM-Newton is significantly improving our understanding of the LMXRB emission processes, accretion geometries and the physical conditions of the material through which the X–ray emission propagates. Reflection Grating Spectrometer (RGS) observations of the eclipsing and dipping LMXRB EXO 0748-676 (Cottam et al. c:01 (2001)) revealed several discrete absorption and emission features from ionized neon (Ne x Ly$`\alpha `$, Ne ix Ly$`\alpha `$), oxygen (O viii Ly$`\alpha `$, O vii He-like complex) and nitrogen (N vii Ly$`\alpha `$). These line features are broadened (velocity widths of 1000–3000 km s-1), but with no evidence for a velocity shift ($`<`$300 km s-1). No obvious dependence on orbital phase was found.

During an XMM-Newton observation of the eclipsing and dipping LMXRB MXB 1659–298 resonant absorption features identified with O viii 1s-2p, 1s-3p and 1s-4p, Ne x 1s-2p, Fe xxv 1s-2p, and Fe xxvi 1s-2p transitions, together with a broad Fe emission feature at $`6.47_{0.14}^{+0.18}`$ keV were discovered (Sidoli et al. si:01 (2001)). The line widths are unresolved in the RGS and correspond to velocities of $`<`$600 km s-1 for O viii 1s-2p and $`<`$2300 km s-1 for the Fe xxvi K$`\alpha `$ line observed using the European Photon Imaging Camera (EPIC). The absorption features occur through a wide range of orbital phases, so the absorbing material is most likely located above, or below, the accretion disk, in a cylindrical geometry with an axis perpendicular to that of the disk.

Similar absorption features have been reported by Parmar et al. (p:01 (2002)) during an XMM-Newton observation of the LMXRB X 1624$``$490, sometimes called the “Big Dipper”. Here the lines are identified with the resonant K$`\alpha `$ absorption of Fe xxv and Fe xxvi, with Ni xxvii K$`\alpha `$ and Fe xxvi K$`\beta `$ possibly detected. Again, no obvious dependence on orbital phase was observed, except during an absorption dip. Here, we report the discovery of discrete X–ray absorption features from highly ionized Fe and Ca (and possibly Ni) in XMM-Newton EPIC spectra of GX 13$`+`$1.

## 2 Observations

The XMM-Newton Observatory (Jansen et al. j:01 (2001)) includes three 1500 cm2 X–ray telescopes each with an EPIC at the focus. Two of the EPIC imaging spectrometers use MOS CCDs (Turner et al. t:01 (2001)) and one uses a PN CCD (Strüder et al. st:01 (2001)). The region of sky containing GX 13$`+`$1 was observed by XMM-Newton a number of times as part of the calibration programme. Three on-axis observations were performed (see Table 1). Since an analysis of two off-axis observations (2000 April 1 13:34 and April 3 13:34 UTC) give similar spectral results, we discuss only the better calibrated on-axis observations. In order to minimize the effects of pile-up the PN and MOS2 were operated in their Timing Modes. The MOS1 CCD was operated in its normal Full Frame Mode and the data are strongly affected by pile-up (when two photons are detected closely in time and the energy is incorrectly assigned to be the sum of the individual energies; see Strüder et al. st:01 (2001)). The effects of pile-up in MOS1 were minimized by extracting events in an annulus outside the 30′′ radius core of the GX 13$`+`$1 point spread function (PSF). Since there are far fewer MOS events than PN ones, the MOS1 data are only used to verify the overall shape and energy calibration of the PN data.

Raw data products were extracted from the public XMM-Newton archive and then reprocessed using version 5.2 of the Science Analysis Software (SAS), before being further filtered using xmmselect. Only X–ray events corresponding to patterns 0 (single pixel events) were selected. Background spectra were obtained from source free regions taken from the same observations. For the Timing Mode observations, appropriate public response matrices are not yet available. Thus, the standard PN matrix for Full Frame Mode has been adopted. This slightly affects the modeling of the continuum, and possibly the energies of the lines (F. Haberl, priv. communication), but is unlikely to affect estimates of their EWs. We concentrate here on the analysis of the narrow absorption features evident in the spectra and defer a detailed discussion of the continuum properties to a later paper.

## 3 Results

A total of 6 EPIC spectra of GX 13$`+`$1 (3 MOS1 and 3 PN) from the on-axis observations reported in Table 1 were extracted. These were then rebinned to oversample the FWHM of the energy resolution by a factor 3 and to have additionally a minimum of 20 counts per bin to allow use of the $`\chi ^2`$ statistic. In order to account for systematic effects a 2% uncertainty was added quadratically to each spectral bin. The photo-electric absorption cross sections of Morrison & McCammon(m:83 (1983)) are used throughout. All spectral uncertainties are given at 90% confidence, unless indicated otherwise.

The overall 2.0–10 keV continua were investigated using the same model as applied to the ASCA spectrum reported in Ueda et al. (u:01 (2001)) consisting of absorbed blackbody and multicolor disk-blackbody (Mitsuda et al. m:84 (1984)) components. Fits were made to all 6 spectra individually. Examination of the fit residuals shows that, especially in the 3 PN spectra, broad positive residuals are present near 6.7 keV, requiring the addition of a broad Gaussian line in all 3 fits with EWs of $``$100 eV and widths, $`\sigma `$, in the range $``$200–400 eV. The continuum parameters derived from all spectra are similar. For the multicolor disk-blackbody component $`kT`$ is in the range 0.5–0.7 keV while the $`kT`$ of the blackbody is 1.1–1.2 keV. Both components are absorbed by an $`N_\mathrm{H}`$, of $`3.1\mathrm{3.8\hspace{0.17em}10}^{22}`$ atom cm-2. Without the addition of the emission and absorption features discussed below the fit quality is rather poor with a $`\chi _\nu ^2`$ of 3.5–4.6. The addition of emission and absorption lines (modeled by Gaussians) successively decreases the $`\chi _\nu ^2`$ to 2.0–2.4 (see below). The spectral parameters for the multicolor disk-blackbody and blackbody are consistent with those obtained from the ASCA observation by Asai et al. (a:00 (2000)).

Two strong absorption features are observed in all 6 spectra at around 6.7 and 7.0 keV (see Table 2). These features are consistent with being due to resonant K shell absorption from highly ionized iron Fe xxv and Fe xxvi, respectively. The deeper 7.0 keV feature is almost certainly the absorption feature observed from GX 13$`+`$1 by Ueda et al. (u:01 (2001)) using ASCA. The addition of both features is highly significant. Including the 7.0 keV line reduces the $`\chi _\nu ^2`$ by 0.6–0.8 while the addition of the 6.7 keV feature reduces the $`\chi _\nu ^2`$ by 0.3–0.6. Two other absorption features at around 7.8 and 8.2 keV are also detected here for the first time in GX 13$`+`$1 in all the PN spectra, while in the MOS1 spectra they are visible only in one observation (Obs. ID 1001), consistent with the lower effective area of this instrument, compared to the PN. The addition of these features reduces the $`\chi _\nu ^2`$ by 0.15–0.3 and 0.04–0.07 for the 7.8 and 8.2 keV features, respectively. In both the MOS1 and PN spectra an absorption feature at $``$4.1 keV is also evident. This is consistent with being due to resonant K$`\alpha `$ absorption from H-like calcium (Ca xx). All the detected absorption lines and their parameters are listed in Table LABEL:tab:lines. We note that there is a strong interplay between the parameters of the broad emission line and the EWs of the two strong absorption features since these features cut strongly into the emission line. This is most clearly seen in Fig. 1 which shows the best-fit PN spectra together with expanded plots of the spectral regions that include the Fe features.

Line energies measured using the PN Timing Mode spectra were expected to be slightly higher than their “true” values due to uncertainties in the calibration of Timing Mode. However, a comparison of the energies derived from the MOS1 and PN spectra for the two strongest absorption lines (Fe xxvi and Fe xxv K$`\alpha `$) shows that they are consistent at 90% confidence. We therefore present in Table LABEL:tab:lines the PN energies without any adjustment. The identification of the fainter absorption features at 7.8 and 8.2 keV, is uncertain. They may be tentatively identified with Fe xxv (7.88 keV) and Fe xxvi K$`\beta `$ (8.26 keV), or with Ni xxvii (7.80 keV) and Ni xxviii K$`\alpha `$ (8.10 keV), respectively (e.g., Kotani et al. k:00 (2000)).

Close inspection of the residual plots (Fig. 1) reveals the presence of an emission feature near 8.5 keV. A fit to this feature with a Gaussian model gives an energy of $`8.51\pm 0.05`$ keV and an EW in the range 25–35 eV. Since this energy is close to those of known PN background features due to Zn K$`\alpha `$ and Cu K$`\beta `$ (an even stronger line is present in the background spectrum near 8.05 keV from Cu K$`\alpha `$) we investigated whether this feature could originate from the background. In the PN Timing Mode, data from the on-axis CCD are continually clocked across the target region and data are extracted with only one dimension of spatial information. We searched all the publically available XMM-Newton observations but were unable to find any blank field Timing Mode data from which to extract a background spectrum. Since the internal background is spatially very structured, particularly in the K- fluorescent lines of metal components of the camera body, we constructed a background spectrum from the identical region covered by Timing Mode from Imaging Mode data, taking care to apply the appropriate exposure time correction. From this spectrum we find that the observed 8.5 keV feature is a factor 100 more intense than the background features. We therefore conclude that the observed feature is not a background artifact. The energy of the emission feature does not correspond to any prominent emission line from an abundant element. However, it is just below that of the He-like Fe edge at 8.83 keV. We therefore attempted to model this feature as an absorption edge. The inclusion of an edge at 8.83 keV accounts for the excess around 8.5 keV and decreases the $`\chi _\nu ^2`$ to 1.6–1.7. The best-fit parameters for all the discrete features can be found in Table LABEL:tab:lines, where the upper limits to the presence of a H-like Fe edge (at 9.28 keV) are also reported. A narrow emission feature near 6.4 keV was reported by Ueda et al. (u:01 (2001)). Addition of this narrow feature with an energy close to 6.4 keV resulted in a decrease of the $`\chi ^2`$ of 4–7, which is not significant.

We note that the parameters derived for the broad emission feature from the PN and MOS1 are inconsistent. The energy derived from the MOS is lower, while the width of the line is significantly broader. This is not caused by the fitting method (e.g. a local minimum) since the addition of the line with the PN-derived parameters to the MOS spectral model results in a clearly worse $`\chi ^2`$ ($`\mathrm{\Delta }`$$`\chi ^2`$$`>`$100), and the difference between the model and the MOS data is up to 10–15%, i.e. larger than the calibration uncertainty. Since a large fraction of the line emission is removed by the narrow absorption components it is difficult to reliably constrain the emission component. We therefore await an improved EPIC calibration before the parameters of the broad emission line can be accurately determined.

## 4 Discussion

There is a growing evidence that narrow absorption lines from highly ionized material are a common feature of LMXRB spectra. They were first detected from the superluminal jet sources GRO J1655$``$40 (Ueda et al. u:98 (1998); Yamaoka et al. y:01 (2001)) and GRS 1915+105 (Kotani et al. k:00 (2000); Lee et al. l:01 (2001)) with ASCA. Absorption lines due to Fe xxv and Fe xxvi were observed from GRO J1655$``$40 and did not show any obvious dependence of their EWs on orbital phase. Both X–ray and optical observations indicate that GRO J1655$``$40 is most probably a binary system viewed at a high inclination angle of $``$$`70^{}`$ (Kuulkers et al. k:98 (1998); van der Hooft et al. vdh:97 (1997); Greene et al. gr:01 (2001)). ASCA observations of GRS 1915+105 revealed, in addition, absorption features due to Ca xx, Ni xxvii and Ni xxviii. This source has been observed also with the Chandra HETGS which detected absorption edges of Fe, Si, Mg, and S as well as resonant absorption features from Fe xxv and Fe xxvi and possibly Ca xx (Lee et al. l:01 (2001)). The recent discovery of similar narrow iron absorption features in X–ray binaries containing neutron stars (GX 13+1, Ueda et al. u:01 (2001); MXB 1658$``$298, Sidoli et al. si:01 (2001) and X 1624-490, Parmar et al. p:01 (2002)), rules out the possibility that such features are peculiar to superluminal jet sources and related in some way to the jet formation mechanism. Chandra HETGS observations have also convincingly detected narrow absorption features due to Ne ix and Fe xviii from the transient black hole candidate XTE J1650$``$500 (Miller et al. m:02 (2002)). These lines are may be due to absorption in an ionized accretion disk atmosphere or wind.

The XMM-Newton observations reported here reveal a much more complex range of emission and absorption components from GX 13$`+`$1 than even the ASCA analysis of Ueda et al. (u:01 (2001)) had found. In addition to the Fe xxvi K$`\alpha `$ absorption feature additional narrow absorption features which may be identified with Ca xx and Fe xxv K$`\alpha `$ and with Fe xxv and Fe xxvi K$`\beta `$ are evident in the XMM-Newton spectra. In addition a deep edge at 8.83 keV is clearly required. We note that the feature at $`7.61\pm 0.13`$ keV which Ueda et al. (u:01 (2001)) model as an edge is more clearly seen to be a line absorption feature in the XMM-Newton PN spectra (which we identify with Fe xxv K$`\beta `$). It is interesting to speculate as to why a deep ($`\tau 0.2`$) edge from Fe xxv is evident whereas the PN upper-limits to the optical depth for an Fe xxvi edge at 9.28 keV of $`\text{ }<`$0.04 exclude a similar feature. This is in contrast to the corresponding absorption lines, where the EW of the Fe xxvi transitions exceeds that of Fe xxv in all the spectra.

In order to convert the observed EWs of the lines into ion column densities, we assumed the curve of growth analysis performed by Kotani et al. (k:00 (2000)). In this way, the column densities of the lines can be estimated for a range of kinematic temperatures of the absorbing material, where this temperature includes contributions from thermal motions as well as any bulk motions or turbulence. An Fe abundance of $`\mathrm{4.7\hspace{0.17em}10}^5`$ was assumed (Anders & Grevesse a:89 (1989)). The column densities for Fe xxv and Fe xxvi, calculated for a range of assumed temperatures (0.1, 1, 10, 100 and 1000 keV), are shown in Fig. 2 for only one observation, for clarity. Indeed, similar EWs have been measured in all the three observations. A lower kinematic temperature requires a higher ion column density, and hence a higher NH. A constraint on the absorbing Fe column density can be derived from the fact that the plasma should be optically thin to Thomson scattering (NH $`<\mathrm{1.5\hspace{0.17em}10}^{24}`$ atom cm-2), otherwise the absorption lines would be strongly diminished. Assuming a solar Fe abundance, implies an Fe column density $`<`$$`\mathrm{7\hspace{0.17em}10}^{19}`$ cm-2, which translates into a temperature $`\text{ }>`$10 keV. The upper limits to the velocity widths of the two narrow lines of $`<`$1800 km s-1 (Fe xxv) and $`<`$1900 km s-1 (Fe xxvi), provide not very constraining upper limits to the $`kT`$ of 600–650 keV. The EWs of the highly ionized Fe K lines (and thus the ions column densities) are similar to those measured in the eclipsing, dipping LMXRB MXB 1658$``$298 (Sidoli et al. si:01 (2001)) with XMM-Newton, and a factor $``$3 higher than that observed in the LMXRB X 1624–490 (Parmar et al. p:01 (2002)). The results of the curve of growth analysis for the Ca xx line are shown in Fig. 3. The ionized calcium column density ranges from 1016 cm-2 to 1020 cm-2, depending on the kinematic temperature. No useful constraints come from the upper limit to the velocity width of $`<`$3700 km s-1, that translates into a temperature of $`\text{ }<`$1900 keV.

The line centroids of the Fe xxv and Fe xxvi K$`\alpha `$ do not show evidence for any velocity shifts, except in the second observation (Obs. ID 0901) where blue-shifts of $`2700\pm 2200`$ km s-1 (Fe xxv K$`\alpha `$) and $`3900\pm 1700`$ km s-1 (Fe xxvi K$`\alpha `$) are evident. In this same observation, a blue-shift in the Ca xx feature of $`5100\pm 2200`$ km s-1 is also detected. For most observations, the MOS energy calibration is correct to within $``$5 eV, or a velocity shift of $``$220 km s-1 at 7 keV and $``$400 km s-1 at 4 keV. Individual observations have been analysed where the Mn K$`\alpha `$ internal calibration line appears to be discrepant by $``$30 eV, but these are all associated with anomalous instrument temperatures resulting from spacecraft irregularities. Such conditions did not pertain for the observations reported here. The absorption lines at higher energy cannot be unambiguously identified. These features are only clearly detected in the MOS1 in the third observation (Obs. ID. 1001). The measured line energy of 7.84$`\pm 0.04`$ keV can be associated with both the Fe xxv K$`\beta `$ (at 7.88 keV) and the Ni xxvii K$`\alpha `$ (at 7.80 keV) feature. The absorption line at $``$8.1 keV can be similarly associated with either an Fe xxvi K$`\beta `$ transition (at 8.212 keV) or with a Ni xxviii K$`\alpha `$ feature at 8.034 keV.

Since the orbital period of GX 13$`+`$1 has not been confirmed, we cannot reliably estimate a geometry and size for the absorbing plasma around the central X–ray source. If the 25 day period is confirmed to be orbital, we note that the XMM-Newton observations reported here only cover a small fraction ($``$10%) of an orbital cycle. All we can say is that all the narrow absorption lines are consistent with being present with the same EW in all analyzed spectra. This implies that the absorbing plasma is stable in the line of sight over timescales of few days suggesting a stable geometry, as in the case of MXB 1658$``$298 and X 1624–490, and excluding the possibility of absorbing clouds temporarily passing in front of the central source producing the absorption. Additionally, the fact that the EWs of the Fe xxvi line measured here are consistent with that observed with ASCA (Ueda et al. u:01 (2001)) a few years ago, again supports this conclusion.

###### Acknowledgements.

Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA member states and the USA (NASA). We thank T. Kotani for making his curve of growth software available.