text:
# FLY. A parallel tree N-body code for cosmological simulations. Reference Guide

## 1 Introduction

FLY is a parallel collisionless N-body code which relies on the hierarchical oct-tree domain decomposition introduced in 1986Natur.324..446B  for the calculation of the gravitational force. Although there exist different publicly available parallel treecodes, FLY differs from them because it heavily relies on two parallel programming concepts: shared memory and one-sided communications. Both of these concepts are implemented in the SHMEM library of the UNICOS operating system on the CRAY T3E and Sgi Origin computing systems.
FLY is the result of the development of a preliminary software called WD99. This code was developed in the period 1996-2000 using several platforms. The first release was developed with the CRAFT programming environment embedded the Cray T3D. The porting of the code on the Cray T3E system, where the CRAFT was no more available, was performed in 1998 using the shmem library, that was the only one-side communication system available. The performances of the shmem library, in terms of scalability and latency time are very good, being this library designed for the hardware architecture of the Cray T3E and of the Sgi Origin systems. On systems like the IBM SP where these libraries are not available FLY has been modified to use the local libraries.
The MPI-2 library was made available with good declared performances for IBM SP System only recently. This library probably will be adopted in the next FLY version, to increase the code portability.
A more detailed treatment of the parallel computing techniques which have been adopted and of the resulting performances on different systems can be found elsewhere.
Being an open source project, FLY can and must be modified to suit the particular needs of individual users. In order to ease this, we will describe some features concerning the integration of the equations of motion (section 2), the generation of initial conditions, the structure of the checkpoint and the output files. In section 5 we will also describe the structure of some parameter files needed to run a simulation. More detailed information concerning the preparation of these parameter files, the compilation and running of the code can be found in the FLY User Guide flyug .

## 2 Equations of motion.

The actual form of the discretized equations of motion is implemented in the subroutines upd\_pos and upd\_vel; if the user needs to identify some special parameters which are used in these subroutines, he may want to define them in initpars. In the publicly distributed version of FLY , we have implemented a set of cosmological equations of motion, solving the standard particle equations of motion for a Friedmann cosmology. In the following, we will present these equations and their discrete implementation.
The Friedmann-Robertson-Walker metric is characterized by an expansion factor $`a(t)`$, where $`t`$ is the conformal time. Let $`𝒙_i(t)`$ be the comoving coordinate of the $`i`$-th particle and $`m_i`$ its mass, then the equations of motion are given by:

$$\dot{𝒙_i}=𝒗_i$$
(1)
$$\dot{𝒗_i}+2\frac{\dot{a}}{a}𝒗_i=\frac{G}{a^3}\underset{ji}{}\frac{m_j(𝒙_i𝒙_j)}{𝒙_i𝒙_j^3}+𝑭_{Ewald}(𝒙)$$
(2)
Here and in the following a dot denotes derivation w.r.t. the conformal time variable $`t`$, $`G`$ is the gravitational constant and the last term, $`𝑭_{Ewald}(𝒙)`$, is the Ewald correction, which takes into account the contribution to the force from the infinite replicas of the simulation box over the spatial directions. We also define the Hubble constant: $`H(t)=\dot{a}/a`$.
It is more convenient to introduce a set of dimensionless spatial, temporal and mass variables:

$$𝒙_i^{}=L_0𝒙_i,t=t_0\tau ,m_i^{}=M_0m_i$$

In terms of these variables, the dimensionless equations of motion become:

$$\frac{d𝒙_i^{}}{d\tau }=𝒗_i^{}$$
(3)
$$\frac{d𝒗_i^{}}{d\tau }+2H(\tau )𝒗_i^{}=\frac{GM_0t_0^2}{L_0^3}\frac{𝒂_i^{}}{a(\tau )^3}$$
(4)
The dimensionless comoving peculiar velocity $`𝒗_i^{}`$ is simply related to the actual peculiar velocity by a scaling relation: $`𝒗_i=(L_0/t_0)𝒗_i^{}`$ (trivially deduced from eq. 3, and the Hubble constant $`H(\tau )`$ in dimensionless time units is given by:

$$H(\tau )=t_0H(t)$$

In fact, the latter relationship is easily obtained after having performed a change of time variable in the definition of the Hubble constant:

$$H(t)=\frac{1}{a}\frac{da}{dt}=\frac{1}{a}\frac{da}{d\tau }\frac{d\tau }{dt}=\frac{1}{at_0}\frac{da}{d\tau }\frac{H(\tau )}{t_0}$$

From now on we will omit the apex from the dimensionless quantities.

### 2.1 Units

We choose units such that:

$$\frac{GM_0t_0^2}{L_0^3}=1$$
(5)
so that the r.h.s. of eq. 4 has a unit multiplying factor. Apart from the obvious reasons of simplicity, this choice helps in reducing the number of floating-point operations in the actual numerical implementation.
We have still to make choices for two out of the three the units entering eq. 5 ($`L_0`$, $`t_0`$ or $`M_0`$). We will adopt the values given by PDBook  for the fundamental constants, rounded to the 6th decimal digit (but in the calculations we use up to the 10th digit when available).
First, note that:

$$GM_{}=1.327125\times 10^{11}\frac{\mathrm{Km}^3}{\mathrm{sec}^2}$$

So, if we measure mass in units of $`M_{}`$ we can rewrite: $`GM_0=GM_{}(M_0/M_{})`$. From now on we will simply write $`M_0`$ in place of $`M_0/M_{}`$.
If we adopt as units of length:

$$L_0=1\mathrm{Mpc}h^1=3.08568\times 10^{19}h^1\mathrm{Km}$$
(6)
and of time:

$$t_0=2/3H_0=2.05759\times 10^{17}h^1\mathrm{sec}$$
(7)
we get from eq. 5 the units of mass:

$$M_0=\frac{L_0^3}{t_0^2GM_{}}=5.22904\times 10^{12}h^1M_{}$$
(8)
These are the fundamental units which we have implicitly adopted in our discretized equations. We can now easily deduce the units of other relevant quantities:

* Velocity:

$$\mathrm{v}_0=\frac{L_0}{t_0}=146.21794\frac{\mathrm{Km}}{\mathrm{sec}}$$
(9)
* Potential (per unit mass):

$$\varphi _0=\left(\frac{L_0}{t_0}\right)^2=2.137971\times 10^4\left(\frac{\mathrm{Km}}{\mathrm{sec}}\right)^2$$
(10)

Finally, we can deduce the values of some useful quantities in these units:

* Gravitational constant. - Using eqs. 5-8, we obtain:

$$\mathrm{G}=\frac{L_0^3}{M_0t_0^2}=\frac{v_0^2L_0}{M_0}$$
(11)
* Hubble constant. - Using the units for velocity we get:

$$\mathrm{H}_0=100\frac{\mathrm{Km}}{\mathrm{sec}\mathrm{Mpc}}=\frac{100}{146.218}\frac{v_0}{L_0}=0.68391\frac{v_0}{L_0}$$
(12)
* Critical Density. - Using the definitions above we get:

$$\rho _c=\frac{3\mathrm{H}_0^2}{8\pi \mathrm{G}}h^2=5.583\times 10^2h^2$$
(13)

### 2.2 Choice of the time variable

In eqs. 3-4 the dimensionless time $`\tau `$ is used as time variable. In actual calculations this is not necessarily the best choice, because often the actual time variable adopted in the outputs is the redshift, which bears a nonlinear relationship with the conformal time, which is also dependent on the underlying cosmological model. We have then adopted a new time variable, introduced first in 1985ApJS…57..241E :

$$p=a^\alpha $$

where $`\alpha `$ is a coefficient to be determined by the user. For closed models ($`\mathrm{\Omega }_0=1`$) G. Efstathiou 1985ApJS…57..241E  advise a value $`\alpha =1/2`$, which ensures an approximately constant r.m.s. advancement of particle positions from one time-step to the next. We choose a constant step $`\mathrm{\Delta }p`$ in order to secure that the leapfrog intergrator keeps second-order accuracy throughout the run.
With this change of variable the new equations of motion are:

$$\frac{d𝒙_i}{dp}=𝒖_i$$
(14)
$$\frac{d𝒖_i}{dp}+2A(p)𝒖_i=B(p)𝒂_i$$
(15)
where the coefficients are given by:

$$A(p)=\frac{1+\alpha +\ddot{a}a/\dot{a}^2}{2\alpha a^\alpha }$$
(16)
$$B(p)=\frac{1}{\alpha ^2\dot{a}^2a^{2\alpha +1}}$$
(17)
(cf. eqs. 10a-b of G. Efstathiou 1985ApJS…57..241E ). The numerical evaluation of the time derivatives $`\dot{a},\ddot{a}`$ is rather cumbersome, but it can be simplified by using some combinations of the Friedmann equations. From S. Carrol (2000, eq. 32) 2000astro-ph…0004075  we have:

$$\ddot{a}a/\dot{a}^2=\mathrm{\Omega }_\mathrm{\Lambda }\frac{1}{2}\mathrm{\Omega }_m$$
(18)
Moreover, one can easily verify (see e.g. 1997astro-ph…9712217 ) that one of the Friedmann equations can be written as:

$$\dot{a}^2a=H^2\left[\mathrm{\Omega }_m(z=0)+\mathrm{\Omega }_ka+\mathrm{\Omega }_\mathrm{\Lambda }(z=0)a^3\right]$$
(19)
Here we have introduced the standard definitions for the dimensionless baryonic ($`\mathrm{\Omega }_m`$), curvature ($`\mathrm{\Omega }_k`$) and vacuum ($`\mathrm{\Omega }_\mathrm{\Lambda }`$) parameters:

$$\mathrm{\Omega }_m=\frac{8\pi \mathrm{G}\rho }{3H^2},\mathrm{\Omega }_\mathrm{\Lambda }=\frac{\mathrm{\Lambda }}{3H^2}$$

Substituting eqs. 18-19 in eqs. 16 and 17, respectively, we get:

$$A(p)=\frac{1+\alpha +\mathrm{\Omega }_\mathrm{\Lambda }\frac{1}{2}\mathrm{\Omega }_m}{2\alpha p}$$
(20)
$$B(p)=\frac{1}{\alpha ^2H^2p^2\left[\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }a^3\right]}$$
(21)
Note that in the equations above the cosmological parameters $`\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda }`$ are all evaluated at $`z=0`$. The numerical evaluation of eqs. 20 and 21 is simpler than that of the original definitions. Some of the factors appearing in these terms are in fact constant and are computed only once in the initialization subroutine inparams.

### 2.3 Remark about initial velocities

As is clear from the preceding paragraph, the velocity $`𝒖`$ does not coincide with the dimensionless peculiar velocity $`𝒗^{}`$. One must keep this fact in mind when computing the initial positions and velocities. From the definition of peculiar velocity in dimensionless units (eq. 3), using the parameter $`p`$ given above, and also eq. 14 we get:

$$𝒗_i^{}=\frac{d𝒙_i}{d\tau }=\frac{d𝒙_i}{dp}\frac{dp}{dt}t_0=𝒖\alpha \left[a(t)\right]^{\alpha 1}\dot{a}\frac{2}{3H_0}=$$

$$=𝒖\alpha \frac{p}{a}\frac{2}{3}a\left[\frac{\mathrm{\Omega }_m}{a^3}+\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/2}$$

Finally, in terms of the expansion factor at the start of the simulation $`a_{in}=(1+z_{in})^1`$, we obtain:

$$𝒖_i=\frac{3}{2\alpha }a_i^\alpha \left[\frac{\mathrm{\Omega }_m}{a_i^3}+\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/2}\frac{𝒗_i}{v_0}$$
(22)
where $`𝒗_i`$ is the peculiar velocity in dimensional units as given for instance, by the COSMICS code 1995astro-ph…9506070

### 2.4 Gravitational potential

The gravitational potential adopted to compute the acceleration is given by the Plummer form

$$\mathrm{\Phi }(r)=\frac{\mathrm{G}m_1m_2}{(r^2+ϵ^2)^{1/2}}$$
(23)
where $`ϵ`$ is a softening length. For $`ϵ>0`$ this potential is finite at the origin, and is adopted in order to avoid the formation of tight binary pairs.

### 2.5 Discretized equations

It is customary to adopt a leapfrog discretization scheme for N-body cosmological simulations 1985ApJS…57..241E , so that our final equations of motion are:

$$𝒖_i^{n+1/2}=𝒖_i^{n1/2}\frac{1A(p_n)\mathrm{\Delta }p}{1+A(p_n)\mathrm{\Delta }p}+\frac{B(p_n)\mathrm{\Delta }p}{1+A(p_n)\mathrm{\Delta }p}𝒂_i^n$$
(24)
$$𝒙_i^{n+1}=𝒙_i^n+𝒖_i^{n+1/2}\mathrm{\Delta }p$$
(25)
where $`𝒂_i^n`$ is the acceleration, the lower index refers to the particle ($`1iN_{bodies}`$) and the upper one to the time-step. Note that positions and velocities are evaluated at steps differing by one-half time-step: this must be taken into account when preparating the initial conditions. Actually before starting the run, the subroutine leap\_corr executes this phase adjustment by advancing the position of half time-step, using a simple second order extrapolation of particles’ positions.

### 2.6 Choice of the time-step

Also the choice of the time stepping criterion should be left to the user’s own choice. In FLY the time stepping criterion is performed in the subroutine dt\_comp. As a default, we adopt a criterion based on the evaluation of the maximum acceleration among all the particles:

$$\mathrm{\Delta }p=\left[\frac{\alpha ϵ_{soft}}{\mathrm{max}_{i=1,N}𝒂_i}\right]$$
(26)
where $`ϵ_{soft},\mathrm{max}_{i=1,N}𝒂_i`$ are the softening length and the maximum acceleration, and $`\alpha 1`$ is a parameter which can be freely adjusted. This time stepping criterion is dynamic, and guarantees an almost constant r.m.s. advancement of the particle with a large velocity. However, the final choice is obviously left to the user, who can easily modify the subroutine dt\_comp or can give a constant time-step (see par. 5.4).

## 3 The FLY grouping

A tree code computes the force on a particle, by means of the interaction between the particle, and some elements (body or cells) of the tree that forms the interaction list $`IL_p`$ for each particle. In a classical tree schema each particle must have its own $`IL_p`$ to compute the force on it.
The fundamental idea of the tree codes consists in the approximation of the force component for a particle. Considering a region $`\gamma `$, the force component on an i-th particle may be computed as

$$\underset{j\gamma }{}\frac{Gm_j𝐝_{ij}}{d_{ij}^3}\frac{GM𝐝_{i,cm}}{d_{i,cm}^3}+\text{higher order multipoles terms }$$
(27)
where $`M=_{j\gamma }m_j`$ and $`cm`$ is the center of mass of $`\gamma `$.
In eq. (27) the multipole expansion is carried out up to the quadrupole order when a far group is considered. The tree method, having no geometrical constraint, adapts dynamically the tree structure to the particles distribution and to the clusters, without loss of accuracy. This method scales as $`O(NlogN)`$.
The criterion for determining whether a cell must be included in the $`IL_p`$ is based on an opening angle parameter $`\theta `$ given by

$$\frac{C_l}{d}\theta ,$$
(28)
where $`C_l`$ is the size of the cell and $`d`$ is the distance of $`p`$ from the center of mass of the cell. The smaller the values of $`\theta `$, the more the number of cells to be opened and hence the more accuracy of forces (for $`\theta =0.8`$ we have an rms error lower than 1% on the accelerations using a multipole expansion of the quadrupole order her87 ).
FLY uses a grouping method to form a grouped interaction list $`IL_g`$, assigned to all the particles lying in grouping cell $`C_{group}`$.
The same above mentioned method is applied to a hypothetical particle placed in the center of mass of the $`C_{group}`$, hereafter VB (Virtual Body). Moreover, we consider the $`IL_g`$ as formed by two parts given by

$$IL_g=IL_{far}+IL_{near}$$
(29)
$`IL_{far}`$ and $`IL_{near}`$ being two subsets of the interaction list. An element is included in one of the two subsets, using the following Sphere method for all the elements that satisfy eq. (28).

Define $`Sphere_{radius}=3\frac{Cellsize(C_{group})\sqrt{3}}{2}`$
If $`Distance(IL_g(element),VB)>Sphere_{radius}`$
Add element to $`IL_{far}`$
Else
Add element to $`IL_{near}`$
Endif

Moreover all the particles $`pC_{group}`$ are included in $`IL_{near}`$.
Using the two subsets it is possible to compute the force $`𝐅_p`$ on a particle $`pC_{group}`$ as the sum of two components,

$$𝐅_p=𝐅_{far}+𝐅_{near},$$
(30)
where $`𝐅_{far}`$ is a force component due to the elements listed in $`IL_{far}`$ and $`𝐅_{near}`$ is the force component due to the elements in $`IL_{near}`$. We assume the component $`𝐅_{far}`$ to be the same for each particle $`pC_{group}`$ and compute it considering the gravitational interaction between the VB and only the elements listed in $`IL_{far}`$, while the $`𝐅_{near}`$ component is computed separately for each $`p`$ particle by the direct interaction with the elements listed in $`IL_{near}`$.

### 3.1 Error Analysis

The dominant component of the error of the classical Barnes-Hut method arise from the truncated multipole series of degree p (generally $`p=2`$) in eq 27, and is O($`\gamma ^{p+1}`$) where here $`\gamma `$ is the inverse of $`\theta `$ (eq. 28) and is always greather than 1. The rms error, in a typical Large scale Structure (LSS) simulation with $`\theta =0.8`$ is about equal to $`1\%`$ on the acceleration, in comparison with a direct particle-particle method 1989ApJS…70..389B .
The error we introduce with the FLY grouping method depends on the choice of some parameters. The ”LIV. GROU” parameter (see par. 5.5) is the level of the tree where a group can be formed. This parameter sets the maximum $`C_{group}`$ size that must be chosen in order to ensure that the difference between $`IL_p`$ and $`IL_{VB}`$ is negligible (no more than $`1\%`$ of the elements). The ”BODY GROU” parameter (see par. 5.5) sets the maximum number of particles that can form a group within a cell: all these particles are listed in $`IL_{near}`$ and there is a direct particle-particle interaction among these elements that form a group.
For a LSS simulation in a 50 Mpc box with more than 2 million particles the ”LIV. GROU” parameter can be set to the sixth level of the tree and the ”BODY GROU” parameter equal to 32. These values maintain the global error around to $`1\%`$ and good performances of the code. Further detail on this method, the errors and the performance of FLY can be found in bec2001 .

## 4 Program organisation

This section list all the main subroutines and functions of FLY . The code is organized in some subroutines that are executed at the start of the job only and other subroutines that are executed for each time-step.

### 4.1 Main and module

fly\_h is the module file of FLY , where all the global variables are defined.
fly is the main program of FLY , that starts the system inizialization and the execution of all the programmed time-step.

### 4.2 System inizialization

The following subroutines are executed before the first time-step, and are listed in the same order of execution during the start-up phase.
null initializes some variables.
sys\_init calls all the following subroutines to read parameters and initial condition files, at the start of the run.
read\_params reads the input parameters from stat\_pars and dyn\_pars files, at the start of the run.
read\_redsh reads the redshift list that will be used to produce the output files.
init\_ew reads and sets the Ewald tables for periodical boundary conditions.
read\_b\_bin reads the binary checkpoint file (i.e. initial condition file) with positions and velocities of all particles, at the start of the run.
read\_b\_asc reads the ascii checkpoint file (i.e. initial condition file) with positions and velocities of all particles, at the start of the run.
init\_balance distributes the initial load among the processors.
init\_pars initializes physical variables at the start of the run.
init\_pos executes the Leapfrog correction on the input data, at the start of the simulation.
reset\_pos resets all particles inside the box where the simulation will evolves, at the start of the simulation.

### 4.3 Time-step execution

The following subroutines are executed for all the time-step of a single job, at the end of the inizialization phase. The number of steps executed by a job being the ”NUM. STEP” value in the stat\_pars file (see par. 5.4). The subroutines are listed in the same order of execution.
inpar\_dyn reads the input parameters from dyn\_pars files, at the start of each time-step cycle.
step executes a time-step cycle: it calls step\_force, dt\_comp, upd\_vel, upd\_pos, out\_32 and wr\_native.
step\_force builds the tree and computes the acceleration on the particles: it calls tree\_build, acc\_comp.
ch\_all allocates RAM for temporary remote data storage.
tree\_build builds the tree data structure. It calls tree\_gen, find\_group and cell\_prop.
tree\_gen builds the geometry of the tree.
find\_group marks all the grouping cell (see sect. 3).
cell\_prop computes the monopole and the quadrupole momentum for all the tree cells.
acc\_comp starts ilist\_group, force\_group, ilist and force subroutines.
ilist\_group walks in the tree and forms the interaction list for each grouping cell.
force\_group computes the grouping acceleration for all particles belonging to a grouping cell, using the interaction list formed by ilist\_group. These components will be completed in the force subroutine.
ilist walks in the tree and forms the interaction list, for each particle.
force computes the acceleration on a particle using the interaction list formed by ilist, or complete the acceleration components for the particles belonging to a grouping cells.
dt\_comp computes the adaptive time-step (see par. 2.6).
upd\_vel advances the particle velocities at the end of a time-step cycle.
upd\_pos advances the particle positions at the end of a time-step cycle.
out\_32 produces the programmed output and the quick-look files using the table read in read\_redsh.
leapf\_corr executes the Leapfrog correction on data output.
wr\_native writes bodies on the checkpoint file: it calls write\_b\_bin or write\_b\_asc.
write\_b\_bin writes the binary checkpoint file, with positions and velocities of all particles, at the end of the run.
write\_b\_asc writes the ascii checkpoint file, with positions and velocities of all particles, at the end of the run.

At the end of each time-step the main program automatically computes the new parameters for the Dynamic Load Balance and writes in the default output the timing of each processor during the step.

## 5 Data input

### 5.1 Initial Conditions

FLY does not give the program for the initial condition generation, but it is possible to use any program the user likes (e.g. COSMICS by E. Bertschinger), provided that the input file is prepared as described here. The input file simply contains all the positions and velocities of all the particles in double precision C format (8 byte for each component). FLY uses a C routine to read and write the input/output particles data file. The initial file input must be produced using the C language, like the FLY io.c function:

```
void write_b_bin(double *pos,double *vel,char *cfilename,signed int *nlong)
{
  FILE *f_b_bin;

  f_b_bin = fopen(cfilename,"wb");
  fwrite(pos,sizeof(double),*nlong,f_b_bin);
  fwrite(vel,sizeof(double),*nlong,f_b_bin);
  fclose(f_b_bin);
}

```

Note that the pos and vel arrays should be defined as: pos\[nbodies\], vel\[nbodies\] in order to keep the same memory map distribution they have in Fortran 90, considering that in the C language, the array ordering is the inverse as in Fortran. So the above ordering corresponds to the following Fortran ordering: pos(3,1:nbodies), vel(3,1:nbodies).
It can often be advisable to create a /temp/FLY directory (where the system allows a great capability area on the /temp file system) to put the initial condition file.
On Cray T3E and Sgi Origin if you use a Fortran program to produce the initial condition file, it should be produced without control characters, i.e. it must have an unblocked file structure with ”stdio” style buffering, compatible with the C fwrite and fread functions. The following code excerpt gives an example:

```
ΨREAL(KIND=8),ΨDIMENSION(3,1:nbodies) :: pos,vel

Ψ....generate the  pos and vel arrays...
Ψ
ΨOPEN(UNIT=12,FILE=’posvel_0’,FORM=’UNFORMATTED’,STATUS=’NEW’)

Ψwrite(12) pos
Ψwrite(12) vel
ΨΨ
ΨCLOSE(UNIT=12)

```

But in the Cray T3E and the Sgi Origin systems you should give the following command at the system prompt, to generate an unblocked file structure with ”stdio” style buffering:

```
prompt>  assign -s sbin u:12

```

Otherwise an ASCII input file can be supplied (for all platforms) with the following data format:

```
ΨREAL(KIND=8),ΨDIMENSION(3,1:nbodies) :: pos,vel

Ψ....generate the  pos and vel arrays...
Ψ
ΨOPEN(UNIT=12,FILE=’posvel_0’,STATUS=’NEW’)

ΨDO i=1,nbodies
ΨΨwrite(12,1000) pos(1:3,i)
ΨENDDO
 ΨDO i=1,nbodies
ΨΨwrite(12,1000) vel(1:3,i)
ΨENDDO

ΨCLOSE(UNIT=12)

1000ΨFORMAT(3(1X,F20.10))

```

### 5.2 Parameter files

You can use the graphical interface to generate all input parameters files, if in the system you use has the Tcl/Tk is installed. The commands

```
prompt> cd ./FLY_2.1/src/tcl
prompt> wish fly_2.1.tcl

```

start the graphical interface that will guide the user to create the parameter files and the Makefile, to compile the code and to create the platform-dependent scripts which allow the user to submit a job in the system queue.
We also provide the user with an assistant program to create the parameter files and to compile FLY . The assistant is invoked using the command assistant\_t3e on Cray T3E (assistant\_ori or assistant\_sp3 on Sgi Origin and IBM SP, respectively) that the user will find in the directory ./FLY\_2.1/src:

```
prompt> cd ./FLY_2.1/src
prompt> ./assistant_t3e

```

The assistant program is structured with several sections and, where applicable, it gives a default value. The user must set a Working Directory and an Executable Directory to store the parameter files and the executable programs respectively. The assistant does not create these directories, and the user must create them before launching the assistant.
After having specified these directories, the assistant will create the input parameter files, and the specific module file that will be used during the compilation.

### 5.3 Generating the ./FLY\_2.1/bin/fly\_fnames file

This file contains the filenames (including the path) of all the input parameter files that will be used during the run. If this file does not exist, the assistant program will generate it and ask for the filenames: the default filenames are reported in square brackets:

```
Starting generation of ../bin/fly_fnames file...

Do you want to create this file (Y/N) : y

All NOT-ABSOLUTE filename paths will be ../bin/ :

Please supply stat_pars filename [stat_pars] :

Please supply dyn_pars filename [dyn_pars] :

Please supply ew_grid filename [ew_grid] :

Please supply ew_tab filename [ew_tab] :

Please supply out32.tab filename [out32.tab] :

Please supply ql.tab filename [ql.tab] :

```

The file will be created, containing the filenames the user gives. In the following we will suppose the user has chosen the default filenames. A more detailed discussion about all the files and the parameters can be found in the FLY User Guide.

### 5.4 The stat\_pars file

This file has the following structure:

```
HEADER   =try_sim
NUM. STEP=  10
MAX_TIME =  0
BAL. PAR.= 0.90
DELTA T. =   0.00450
DT VAR.  =T
OPEN PAR.= 0.90
SOFT PAR.=   0.01000
QUADRUP. =T
OMEGA.CDM=   0.30000
OMEGA.HDM=   0.30000
OMEGA.LAM=   0.40000
HUB_CONST=  0.650000
N. BODIES=     2097152
MASS BODY=   0.00120
IBOD_FILE=/tmp/FLY/posvel_
IBOD_TYPE=B
OBOD_FILE=/tmp/FLY/out_
OBOD_TYPE=B
QLK_FILE =/tmp/FLY/qlk_

```

This file contains some data on the simulations that must be never changed during the system evolution. The assistant asks for each parameter and, eventually, gives a default value. Some peculiar parameters are:

* HEADER is a string to identify the simulation.
* NUM. STEP is the number of time-step cycles of each job.
* MAX\_TIME is the maximum CPU time allowed for each job. This value must be lower than or equal to the single CPU time as specified in the script job. If this value is greater than 0 the NUM. STEP parameter is negligible and the number of time-step cycles for a single job, is automatically computed by FLY .
* BAL. PAR. is used to balance the load among the processors at the start of a run. This value is automatically computed at the end of each time-step cycle, to allow a perfectly load balance. It is the percentage of the local particles that must be computed by the local processor, the remaining portion being computed by the first available processor. We recommend to use the default value 0.90.
* DELTA T. is the integrator time-step used in the system evolution . For cosmological simulations a possible choice can be found in 1985ApJS…57..241E .
* DT VAR. FLY can also use an adaptive Delta T integrator automatically computed in dt\_comp.F subroutine (see par. 2.6). In this case the user must give the T (True) value to this parameter. The F (False) value, forces FLY to use a fixed Delta T as given above.
* OPEN PAR is the opening angle parameter ($`\theta `$), used to open or close the cells during the tree walking phase, to build the interaction list for each particle (eq.28).
* SOFT PAR. is the softening length of the gravitational interaction (see par. 2.4).
* IBOD\_FILE This is the root filename of the checkpoint file, to which a counter is automatically added (as a suffix) to form the complete filename.
* IBOD\_TYPE. The user must specify B for binary or A for ascii data format for the IBOD\_FILE data type.
* OBOD\_FILE This is the root filename of the output files, to which a counter, giving the redshift value, is automatically added to form the complete filename. These output files are the scientific data output produced by FLY with the Leapfrog correction for the phase adjustment. The data format of these files is the same as the checkpoint files.
* OBOD\_TYPE. The user must specify B for binary or A for ascii data format for the OBOD\_FILE data type.
* QLK\_FILE This is the root filename of the quick-look files, to which a counter, giving the redshift value, is automatically added to form the complete filename. This is the ASCII file containing random particles positions and velocity( see par.5.6). During the system simulation a header file (i.e. /tmp/FLY/qlk\_hea) will be automatically created. This file can be used to visualized the system evolution using this file as the header file of AstroMD (http://www.cineca.it/astromd)

The qlk\_hea file has a format like the following:

```
header_line
ASCII
time
262144
17  /gpfs/temp/ube/data262k/qlk_60.0000
10  /gpfs/temp/ube/data262k/qlk_50.0000
68  /gpfs/temp/ube/data262k/qlk_20.0000

```

The value 262144 is the number of data points included in the quick-look file, the filename suffix ( \_60.0000) is the corresponding redshift value when the output is produced, and the number before the filenames (17, 10, 68) is the number of time-step cycles that FLY has reined from the start or the previous quick-look file produced.

### 5.5 The dyn\_pars file

This file contains some parameters of a simulation, that could be changed during the system evolution, and has the following structure:

```
CURR.TIME=  50.000000000000000
CURR.STEP=  0
MAX  STEP=100
LIV. GROU=  9
BODY GROU= 16
GROUP FL.=  1
SORT_LEV.=  3
BOX COMP.=  0
BOX SIZE =  50.0000
X MIN VER=  0.000000
Y MIN VER=  0.000000
Z MIN VER=  0.000000

```

This file must be created by the user and is read at the beginning of each time-step cycle during a run. This file is automatically updated at the end of each job. The assistant asks for each parameter and, eventually, gives a default value. Some peculiar parameter are:

* CURR.TIME is the redshift value at the current time-step cycle.
* CURR.STEP. is the current time-step cycle. We strongly recommend to start from 0 as initial time-step cycle
* MAX STEP. is the maximum allowable time-step cycle. FLY halts the simulation when this value, or the final redshift given by the user (see par. 5.6) is reached.
* LIV. GROU is the level where the grouping cell can be started to form, being LIV. GROU=0 the root level .
* BODY GROU is the maximum number of particles $`N_{crit}`$ that a grouping cell can contain (see above). We recommend to use a value from 8 up to 32.
* SORT LEV. is automatically computed by the assistant, and is used by the FLY\_sort utility, that builds a tree up to the level indicated by this parameter, and produces a sorted data input file, organized as the tree structure

The items GROUP FL. and BOX COMP are for future usage, and the items BOX SIZE, X-Y-Z MIN VER are respectively the size of the box and the coordinates of the lowest vertex.

### 5.6 Other files

The following files must be located in the working directory of FLY and contain the following data.

* out32.tab. List of programmed redshift outputs, FLY generates the output and the quick-look files, with the root filename reported in the OBOD\_FILE and QLK\_FILE parameters of the stat\_pars file, and a suffix given by the redshift value listed in this table.
* ql.tab is a random sequence of values used to generate quick-look file.
* ew\_grid and ew\_tab. These tables are used for the Ewald correction considering the boundary periodical conditions.

More details on these files can be found in the FLY User Guide.

## 6 Installing FLY

The FLY installation procedure is straightforward and even the most unexperienced reader can compile and run successfully a parallel code.
Download the file FLY.tar.gz and copy it somewhere (e.g. in your HOME directory) of a Cray T3E system and/or Sgi Origin and/or an IBM SP system.
Make sure that a subdirectory with the same name ./FLY\_2.1 is not present, to avoid possible conflicts.

### 6.1 Unpacking FLY

Once you have downloaded the package, give the following commands to uncompress and unfold it:
$ gunzip FLY.tar.gz
$ tar -xvf FLY.tar

The tar command will create the following directories:

./FLY\_2.1 (FLY installation directory)
./FLY\_2.1/bin (Executable programs and input parameters files)
./FLY\_2.1/src (Source code and assistant program)
./FLY\_2.1/src/tcl (TCL program: graphical interface)
./FLY\_2.1/job (Job script to run FLY)
./FLY\_2.1/out\_stat (Output log file of the FLY run)
./FLY\_2.1/out\_log (Output log file of the job script)
./FLY\_2.1/Testcase (Input/output example files)
./FLY\_2.1/tools (Utilities)
./FLY\_2.1/FLY\_ug.doc (FLY User Guide)

Now you can use the graphical interface to generate all the parameters files and the executable program (see par. 5.2).

### 6.2 FLY module, makefile and script files

The fly\_h.F file is a module file containing all the fundamental data that are common to all the subroutines of FLY . To generate an executable program, the user must set some fundamental parameters.
The user must give the number of processors of the parallel execution, the number of particles that he wants to use (the same as stat\_pars file), and the estimated ram available for each processor (this value is used to allocate dynamically a temporary buffer in the local ram), and the maximum length of the interaction list $`IL`$, formed during the tree walk procedure: this value depends on the box size, the number of bodies, and the opening parameter: in a clustered regions this value is greater than in an uniform region. It is very important, for the FLY performance, to give a safety value but close to the expected maximum length of the interaction list.
All the parameters have a default value that is set taking a standard cosmological model as a template.
The mkfl\_fly file is the makefile of FLY . The user must give the filename of the executable program and the executable directory, the user can also specify if he wants receive a statistical report. The statistical report concerns only the code performance, it must be used only from expert users of FLY , and must not be used during a production run. FLY also has the mkfl\_fly\_sort file to generate the FLY\_sort utility.
The last file the user can create is the FLY\_job.cmd. This is a script file that contains a schema to submit a job to the system queue. FLY also creates the FLY\_sort\_job.cmd to submit the FLY\_sort utility. On IBM SP, some information about the node characteristics and the system topology are also asked to generate the script file. In any case the user must customize both the FLY\_job.cmd and FLY\_sort\_job.cmd files.

#### 6.2.1 Compiler options

FLY automatically sets the compiler options, that can be changed by expert users. All the options were well tested in order to reach the best performance of the code. The following options are given:

Cray T3E system

* compilers: f90 and cc
* -N132 : specify the line width to 132 column line;
* -O3 : aggressive optimization;
* -LANG:recursive=on : the compiler assumes that a statically allocated local variable could be referenced or modified by a recursive procedure call;
* -DT3E : must be always given. The T3E flag is used to compile specific section of the code for this specific platform;
* -DPW2 : must be used to obtain higher performances when the number of processors that will be used to run a simulation is a power of two;
* -DSTA : is used to obtain a statistical report on the number of remote GET and PUT calls. can be used to test the fly efficiency and must not be used during the simulation run;
* -DSORT : generate the FLY \_sort utility.

SGI Origin system

* compilers: f90 and cc
* -64 : generate a 64 bit objects;
* -extended\_source : specify the line width to 132 column line;
* -O3 : see Cray T3E system option;
* -LANG:recursive=on : see Cray T3E system option;
* -DORIGIN, -DPW2, -DSTA, -DSORT : see Cray T3E system options;

IBM SP system

* compilers: mpxlf90\_r and xlc
* -O4 : aggressive optimization;
* -qrealsize=8 : sets default real to 8 bytes;
* -qfloat=fltint:rsqrt : optimizations for floating point operations;
* -qmaxmem=-1 : unlimited memory used by space intensive optimizations;
* -DSP3, -DPW2, -DSTA, -DSORT : see Cray T3E system options;

## 7 Starting a simulation

The parameters that the user must set to start a cosmological simulation, are listed in the stat\_pars file (see par. 5.4). The most important parameters that describe the cosmological model and the time evolution of a simulation are listed in the following.

* DELTA T. Is the integrator time-step used in the system evolution. FLY uses the Leapfrog second order integrator for the system dynamic. The user can set a constant safe value to ensure the accuracy, the stability and the efficiency of the run. This value is ignored if the user want to use the adaptive integrator time-step, included in the FLY code.
* DT VAR. This is a logical value. A True value (T) means that the user want use an adaptive integrator time-step. FLY computes each time-step the integrator value as described in the eq. 26
* OPEN PAR. This value is the $`\theta `$ opening angle parameter that is used in the Barnes-Hut algorithm. A discussion of this parameter is given in the sect. 3 eq. 28. Typical values used for the LSS cosmological simulation are in the range 0.6 - 1.0.
* SOFT PAR A softened force of the form

$$𝑭=\frac{GMm𝑹}{(R^2+ϵ^2)^{3/2}}$$
(31)
is adopted in the FLY code, where $`ϵ`$ is the softening parameter.
This parameter is crucial in determining the amount of hardening of the gravitational interaction ultimately the amount of small scale substructure which is not destroyed during the gravitational evolution of the system.
* QUADRUP. This is a logical value. A true value (T) means that the quadrupole order will be used in the eq. 27. A false value means that the higher order multipole terms of the force component will not be considered. In order to have a good accuracy for the force component, the true value is strongly recommended.
* OMEGA CDM, OMEGA HDM and OMEGA LAM These parameters must be set by the scientist and are used to select a cosmological model for the simulation. The summation of all these values must be equal to 1.0.
* HUB\_CONST Is the Hubble constant $`H_0`$ (eq. 12) express in the unit of $`100Km/sec/Mpc`$. A value of 0.65 means $`H_0=65Km/sec/Mpc`$.
* N. BODIES Is the total number of particles of a simulation.
* MASS BODY Is the mass of each particle, the unit of this quantity is given in eq. 8.

Other parameters that can be changed during the run, are listed in the dyn\_pars file (see par. 5.5).

* LIV GROU and BODY GROU are described in the par. 3.1.
* BOX SIZE Is the size of the box region where the LSS simulation is executed.
* X - Y - Z MIN VER represents the cartesian coordinates of the lowest vertex of the box. All the particles of the simulation must be included in the box.

The remaining parameters of the stat\_pars and dyn\_pars files do not regard the physics of a simulation, but only the I/O filenames and the job duration, and are already described in the par. 5.4 and 5.5.

## 8 Running FLY

The following procedure is recommended to perform a simulation using several FLY runs:

* Create the input file at the time-step cycle equal to 0 (i.e. /tmp/FLY/in\_cond\_0). We suggest to run the FLY\_sort utility that will generate /tmp/FLY/in\_cond\_sort\_0, the new sorted input file, and replace it as the new input file (i.e. /tmp/FLY/posvel\_0).
* Submit the FLY job using the FLY\_job command file. We suggest to use the ckstop file that is automatically created by FLY when the executable stops the programmed time-step cycles, but the simulation has not yet reached the final redshift.

The User Guide gives more details about the use of the command files.

## 9 FLY scalability

In this section we report scalability data of the FLY code using a testcase of 2097152 particles in a clustered configuration using the Sgi Origin 3800 system at the Cineca. It is a 128 processor elements (PEs) RISC 14000, 500 Mhz with 2 GBytes DRAM for each processor. The Fig. 1 shows the scalability, using the system from 2 PEs up to 30 PEs, the maximum PE number available for each run on this system. The Fig. 2 shows the scalability of the code. We fix the optimal number of processor to 16 and we increase the number of particles in the simulation, from 1 million up to 32 million-particles. All the runs were carried out using $`\theta =0.8`$, the quadrupole expansion, and a grouping level not lower than 6. More details on the FLY scalability are reported in bec2001 .

## 10 Testcase

In order to check that FLY is working correctly after it has been installed, the user can use the testcase reported in the FLY distribution (Testcase direcory). Change the subdirectory of the system platform where you are running (i.e. ./FLY\_2.1/Testcase/SP) and create the FLY executable in this directory to run a test simulation with 1024 particles and 10 time-step cycles. Use the supplied files: fly\_fnames, stat\_pars, dyn\_pars, ew\_grid, ew\_tab, ql.dat, out32.tab, the input ascii file 1024\_asc\_10 and submit the job. The execution will produce the 1024\_asc\_20 file that must be equal to output\_gr file included in the FLY distribution. The testcase also gives another testcase with 2097152 particles: use the fly\_fnames\_2M (copy it as fly\_fnames) and create the FLY executable to run it with this number of particles. In this case FLY must run without errors and the tree properties, the number of cells formed level-by-level, and the root cell properties (written by the FLY run in the standard output) must be equal to the output\_tree\_2M included in the Testcase directory.

## 11 FLY graphical interface

All the input parameter files can be created using the on-line assistant that does not need a graphical environment. Moreover if the user have installed the wish Tcl/Tk in the system where he want to run FLY , he can use the graphical interface in ./FLY\_2.1/src/tcl. He must run wish fly\_2.1.tcl command to start the graphical interface that help the user to create all the parameter files, excluding the initial condition file.
The first time that this interface is started, it asks for the platform, this data and other default data will be saved in the fly\_tcl.ini file. The main window sets the working directory, the executable directory and will create the not existing directories. The figures 1 and 2 show the main window and the window to set the stat\_pars file. If this file exists, it is loaded with the values of the existing file. Other files are created using similar windows.

In the main window the user can click on the Generate button and insert data to create the fly\_h.F module, the makefiles in the directory src, and the script file that can be submitted to the system queue. The Make button executes the makefile and creates the executable program and the FLY\_sort utility.
In the main window it’s also possible to click on the buttons Interactive Run and Batch Run to submit the FLY execution.

## 12 FLY tools

FLY can produce very large output files, whose analysis could be long and cumbersome even on large computing systems. We have then included two utilities which help the user to make some simple analysis using widespread visualization packages, like SuperMongo and IDL. They are contained in the subdirectory tools, and are essentially a C subroutine, slide\_3\_2d.c, an IDL procedure (slice.pro) and a SuperMongo macro: pl\_qlk.
The program slide\_3\_2d.c is an interactive utility which reads a quick-look input and outputs a file containing densities computed on a grid which is a 2-D projection of a slice of the system along the chosen line-of-sight. The latter is specified as an input parameter. The output is a binary file (slice.d) containing the densities on a grid. For instance, in order to get from a quick-look file qlk\_3.0, containing $`262144`$ positions, a slice along the $`y`$direction, on a grid containing $`10^3`$ cells using only particles having $`15<y<22.4`$ the command would be:

$$\$\mathrm{slices2d}\mathrm{n}\mathrm{qlk}\mathrm{\_}3.0\mathrm{nb}\mathrm{\hspace{0.17em}262144}\mathrm{nc}\mathrm{\hspace{0.17em}10}\mathrm{d}\mathrm{\hspace{0.17em}2}\mathrm{min}\mathrm{\hspace{0.17em}15}$$

$$\mathrm{max}\mathrm{\hspace{0.17em}22.4}$$

The IDL procedure slice.pro makes isodensity contours of the file slice.d, while the SMongo macro pl\_qlk plots the particles on a given slice.
Another tool is the tofly.c routine executing a conversion from big-endian to little-endian and viceversa.

## 13 Troubleshouting

Using the assistant and the graphical interface, the user sets parameters of the fly\_h.F module and uses it to create the FLY executable. But many other parameters must be adjusted by the user if some errors occur. This section describes the main errors that can occur and the action that must be taken to avoid the error. All the following parameters are included in the fly\_h.F module. After any change in the parameter value, FLY must be re-compiled.

* ch\_all error 1: Abort. It is not possible to have the minimum RAM-cache size. Action: there is not enough available memory for this execution. Increase the available memory if it is possible, or decrease the nb\_loc (number of local body) or nc\_loc (number of local cells) parameter value.
* find\_group error 1: overflow. Action: increase the maxilf parameter value, the maximum length of temporary storage to form an Interaction List.
* ilist error 1: overflow. Action: increase the maxnterm parameter value.
* ilist error 2: overflow. Action: increase the maxilf parameter value. If the error will persist, please report the error to the FLY authors.
* ilist\_group error 1: overflow. Action: increase the maxnterm parameter value.
* ilist\_group error 3: overflow. Action: increase the maxilf parameter value. If the error will persist, please report the error to the FLY authors.
* read\_params error 1: nbodies greater than nbodsmax. Action: the number value nbodies read from the stat\_pars file doesn’t match the nbodsmax parameter. Set nbodsmax equal to nbodies.
* tree\_gen error 1: Max level reached: lev greater then nmax\_level. Action: The number of levels created to build the tree structure is greater than the maximum allowable levels (equal to lmax). Increase the value of lmax parameter.
* tree\_gen error 2: overflow. Action: The number of cells created to build the tree structure, is greater than the maximum allowable internal tree cells (equal to nbodsmax). Increase the value of ncells parameter.

A complete description of all the errors that can occur is reported in the FLY User Guide.

## 14 Conclusion

FLY is still under development, so many features will be added in the future version. Any question, problem and bug can be reported to the authors sending a detailed e-mail to fly\_admin@ct.astro.it, giving a standard test problem, the input parameter files and a description of the system where the bug was found (i.e. operating system, platform, available memory, etc.). An User Guide of FLY 2.1 is available with the FLY distribution flyug . The FLY distribution includes a testcase for each supported platform, allowing the user to start a simple demo of a FLY run. The next public version of FLY will include the computation of the gravitational potential in an adaptive mesh like Paramesh 1999AAS…195.4203O  that will allow the user to interface FLY outputs with hydrodynamic codes that use adaptive mesh. Moreover a FAQ will be prepared and will be accessible from the FLY site http://www.ct.astro.it/fly.


format: text

text:
# Out of the Chorus Line : What Makes 1E1207.4-5209 a Unique Object?

## 1. Introduction

Neutron star atmosphere models predicted the presence of absorption features depending on atmospheric composition, but high quality spectra, collected both by Chandra and by XMM-Newton , did not yield evidence for any feature (see Pavlov et al. 2002a and Becker and Aschenbach 2002 for recent reviews). INS spectra are well fitted by one or more black-body curves with, possibly, a power law contribution at higher energies, but with no absorption or emission features.
The spectrum of 1E1207.4-5209, on the contrary, is dominated by two broad absorption features seen, at 0.7 and 1.4 keV, both by Chandra (Sanwal et al, 2002) and XMM-Newton (Mereghetti et al, 2002). To better understand the nature of such features, XMM-Newton devoted two orbits, for a total observing time of 257,303 sec, to 1E1207.4-5209. In the two MOS EPIC cameras the source yielded 74,600 and 76,700 photons in the energy range 0.2 - 3.5 keV, while the pn camera recorded 208,000 photons, time-tagged to allow for timing studies. Analysis of this long observation, while confirming the two phase-dependent absorption lines at 0.7 and 1.4 keV, unveiled a statistically significant third line at $``$2.1 keV, as well as a possible fourth feature at 2.8 keV. The nearly 1:2:3:4 ratio of the line centroids, as well as the phase variation, naturally following the pulsar B-field rotation, strongly suggest that such lines are due to cyclotron resonance scattering (Bignami et al. 2003). A recent software release, based on a better characterization of the EPIC instrument, prompted us to revisit the data. While the spectral analysis results confirm and strengthen the conlusions of Bignami et al. (2003), the temporal and spatial analysis yielded interesting new results which we shall briefly outline (see De Luca et al. 2003 for details).

## 2. Timing Analysis

After converting the arrival times of the 208,000 pn photons to the Solar System Barycenter, we searched the period range from 424.12 to 424.14 ms using both a folding algorithm with 8 phase bins and the Rayleigh test. The best period value and its uncertainty (P = 424.13076$`\pm `$0.00002 ms) were determined following the procedure outlined in Mereghetti et al. (2002). Comparing the new period measurement of 1E 1207.4–5209 with that obtained with Chandra in January 2000 (Pavlov et al. 2002b), we obtain a period derivative $`\dot{P}`$=(1.4$`\pm `$0.3)$`\times `$10-14 s s-1. However, Fig.1 (left panel) shows that the $`\dot{P}`$ value rests totally on the first Chandra period measurement. Using only the 3 most recent values, the period derivative is unconstrained. Thus, we cannot exclude that the observed spin-down, based on only a few sparse measurements, be affected by glitches or Doppler shifts induced by orbital motion. Questioning the object’s $`\dot{P}`$ would have far reaching consequences for the understanding of 1E1207.4-5209 since the serious discrepancy between the pulsar characteritic age ($`\tau _c5\times 10^5`$ yrs) and the SNR age ($`\tau _{SNR}7`$ kyrs) is entirely based on the value inferred from the measurements summarized in Fig. 1 (left panel).

To study the energy dependence of the pulse profile, we divided the data in four channels with approximately 52,000 counts each: 0.2-0.52 keV, 0.52-0.82 keV, 0.82-1.14 keV and 1.14-3.5 keV. The pulse profiles in the different energy ranges (Fig. 1, right panel) show a broad, nearly sinusoidal shape, with a pulsed fraction varying from $``$ 3 to $``$ 11 % in the four energy intervals. It is worth noting that the minimum pulsed fraction is found in the 0.20-0.52 keV energy range, the only portion of the spectrum free from absorption lines. Indeed, Fig. 1 (right panel) is an independent confirmation of the findings of Bignami et al (2003) who ascribed the source pulsation to the absorption lines phase variation.

Finally, comparing the shapes of the light curves of Fig. 1 (right), we see for the first time a phase shift of nearly 90 between the profile in the lowest energy range ($`<`$0.52 keV) and those at higher energies.

## 3. Optimizing the X-ray position

To derive the sky coordinates of 1E 1207.4–5209 we computed independently for the MOS1 and MOS2 cameras the boresight correction to be applied to the default EPIC astrometry. We used the Guide Star Catalog II (GSC-II111http://www-gsss.stsci.edu/gsc/gsc2/GSC2home.htm) to select, amongst our $``$200 serendipitous detections, 6 sources with a stellar counterpart to be used to correct the EPIC astrometry. The rms error between the refined X-ray and GSC-II positions is $``$1 arcsec per coordinate. The resulting MOS1 position of 1E 1207.4–5209 is $`\alpha _{J2000}=12h10m00.91s`$, $`\delta _{J2000}=52^{}26^{}28.8\mathrm{"}`$ with an overall error radius of 1.5 arcsec. The MOS2 position is $`\alpha _{J2000}=12h10m00.84s`$, $`\delta _{J2000}=52^{}26^{}27.6\mathrm{"}`$, with an uncertainty of 1.5 arcsec, fully consistent with the MOS1 coordinates.

## 4. Search for the optical counterpart

The field of 1E 1207$``$5209 was observed with the 8.2-meter UT-1 Telescope (Antu) of the ESO VLT (Paranal Observatory). Observations were performed with the the FOcal Reducer and Spectrograph 1 (FORS1) instrument. Images were acquired through the Bessel $`V`$ and $`R`$ filters for a total integration time of $``$2 and 3 hrs, respectively. Fig. 2 shows the inner portion of the combined FORS1 $`V`$-band image centered on the target position, with the MOS1 and MOS2 error cicles superimposed. A faint object (marked with the two ticks in Fig. 2) is detected just outside the southern edge of the MOS1 error circle and showed variability along the time span covered by our observations. In any case, its position falls more than 2 arcsec away from the intersection of the MOS1/MOS2 error circles, which we regard as the most probable region.
No candidate counterpart is detected within the MOS1/MOS2 error circles down to R$``$27.1 and V$``$27.3, which we assume as upper limits on the optical flux of 1E 1207$``$5209. For the X-ray derived interstellar absorption ($`A_V=0.65`$) and a distance of 2 kpc, the measured u.l. rule out any hypothetical “normal” stellar companion other than a very low-mass main sequence star. If we assume that 1E 1207-5209 is indeed isolated, we can derive a neutron star optical luminosity $`3.410^{28}`$ erg s-1 or $`4.610^6`$ of its rotational energy loss, a value similar to those of middle-aged INSs. Since the VLT flux upper limits are $`100`$ higher than the extrapolation of the XMM-Newton  blackbody (De Luca et al. 2003), they can not constrain the optical spectrum.

## References

Bignami, G.F., Caraveo, P.A., De Luca, A., Mereghetti, S. 2003 Nature, 423, 72

Becker W. and Aschenbach, B. 2002, in Neutron Stars, Pulsars, and Supernova Remnants, W. Becker et al. eds., p.64

De Luca, A., Mereghetti, S., Caraveo, P.A., Moroni, M., Mignani, R.P., Bignami, G.F. 2003, A&A, submitted

Mereghetti, S., De Luca, A., Caraveo, P.A., Becker, W. Mignani,R., Bignami, G.F. 2002, ApJ, 581, 1280

Pavlov, G.G., Zavlin, V.E., Sanwal, D. 2002a, in Neutron Stars, Pulsars, and Supernova Remnants, W. Becker et al. eds., p.273

Pavlov G.G., Zavlin, V.E., Sanwal, D., Truemper, J. 2002b, ApJ, 569, L95

Sanwal D., Pavlov G.G., Zavlin V.E & Teter M.A. 2002, ApJ, 574, L61


text:
# Cessation of X-ray Pulsation of GX 1+4

## 1 Introduction

Accreting X-ray pulsars are neutron stars in a binary configuration. The X-ray emission of such a system is powered by the accretion of matter from the companion star onto the neutron star. Therefore, unlike their isolated counterparts, neutron stars in a binary configuration are strongly influenced by the mass accretion process.

In many cases, it is believed that an accretion disk is formed around the neutron star in an accreting pulsar. The disk is disrupted by the strong magnetic field of the neutron star at the magnetospheric radius, where plasmas in the disk are channeled to the polar regions of the star along the field lines, creating “hot spots” there. However, if the magnetospheric radius is larger than the co-rotation radius of the neutron star, the accreted matter hits the centrifugal barrier and is likely to be expelled from the system (Illarionov & Sunyaev 1975; Wang & Robertson 1985; Lovelace, Romanova, & Bisnovatyi-Kogan 1999; Bogovalov 2001). This phenomenon is often referred to as the propeller effect.

Observationally, the propeller effect may manifest itself directly in the cessation of pulsation of an accreting pulsar, when the mass accretion rate is sufficiently low (and thus the magnetosphere is sufficiently large). This was indeed first observed in two systems, GX 1+4 and GRO J1744-28 (Cui 1997), which possess contrasting magnetic properties. It was seen again recently, based on our monitoring observations of GX 1+4 during its transition to a rare X-ray faint state (which was very similar to the 1996 episode; Cui 1997). In this paper, we present the new data. We show that the results obtained are in good agreement with those of previous measurements and thus provide additional support for the earlier interpretation.

## 2 Observations

We monitored GX 1+4 weekly with the Rossi X-ray Timing Explorer (RXTE) over a two-year period (2001–2002). Typical observing times are 2–3 ks (roughly one RXTE orbit), although some longer observations were taken, in coordination with the Chandra X-ray Observatory, following the triggering of a Target-of-Opportunity program. For this work, we will focus on the RXTE observations taken in 2002, when a transition to a faint state occurred. The results of the Chandra observation will be discussed elsewhere.

There are two detectors aboard RXTE (Bradt, Rothschild, & Swank 1993), the Proportional Counter Array (PCA) and the High-Energy X-ray Timing Experiment (HEXTE). The PCA consists of five nearly identical large area proportional counter units (PCUs). It has a total collecting area of about 6500 $`cm^2`$ and covers a nominal energy range of 2–60 keV. The HEXTE has a total effective area of about 1600 $`cm^2`$ in two clusters. Each cluster contains 4 detectors, although one of the detectors lost its spectral capability very early in the mission. The HEXTE covers a wide energy range of about 15–250 keV, overlapping the PCA passband at low energies. Since the HEXTE clusters alternate between on-source and off-source fields, the effective exposure time (with dead-time subtracted) is always much less than the actual observing time.

Operational constraints require that some PCUs be turned off. Which PCUs are turned off varies from observation to observation. In our case, the number of PCUs that were operational varied from 2 to 5, with PCU 0 and PCU 2 being on throughout the entire campaign. Because PCU 0 has lost its front veto layer, which makes background estimation more uncertain, we chose PCU 2 as our “standard” detector to normalize the measured fluxes.

## 3 Data Reduction and Analysis

We used the software package HEASOFT 5.2,111see http://heasarc.gsfc.nasa.gov/docs/software/lheasoft along with the associated calibration files and background models, to reduce the RXTE data for spectral analysis. Specifically, Ftools 5.2 was used to extract spectra, both from PCA data and HEXTE data, and XSPEC 11.2 was used for spectral modeling. For timing analysis, however, we used custom software developed at MIT. The software is mature and has been used by the ASM team since the beginning of the mission. It was designed to work directly with raw packet data, instead of standard FITS files.

### 3.1 X-ray Spectral Analysis

To minimize calibration uncertainties associated with the PCA response matrices, we limited ourselves to data from the first xenon layer of each PCU that was turned on during an observation. The trade-off is that we lost sensitivity to X-rays at energies roughly above 25 keV. We used the HEXTE data to constrain the higher-energy portion of the X-ray spectra. It should be noted that cross flux-calibration between the PCA and the HEXTE remains an unresolved issue.222see http://lheawww.gsfc.nasa.gov/docs/xray/xte/crosscal. The effect was taken into account in joint PCA/HEXTE spectral fits by making the relative normalization of the detectors an additional free parameter (i.e., multiplying the model by a constant in XSPEC). This approach has been commonly adopted (e.g., Coburn et al. 2002). Moreover, we chose to extract a spectrum from each PCU and each HEXTE cluster separately and to perform a global fit for each observation, while floating the relative normalization of each spectrum (with respect to PCU 2). This allowed us to also take into account more subtle differences in the effective areas of individual detectors. It also added flexibility in our dealing with the PCUs being turned on and off in the middle of an observation.

### 3.2 X-ray Timing Analysis

To derive properties of the pulsed emission, we made a light curve for each observation from the Standard 1 data, which has a time resolution of 1/8 s and covers the entire PCA passing band. The light curve was then barycenter-corrected and folded at trial periods close to the expected value (since the pulse period is not expected to vary significantly from week to week). At each trial period ($`P`$), we compute the following statistic (Leahy et al. 1983):

$$S=\underset{i=1}{\overset{N}{}}\frac{(r_i\overline{r})^2}{\overline{r}}t_i,$$
(1)
where $`N`$ is the number of phase bins, $`r_i`$ is the average count rate in phase bin $`i`$, $`t_i`$ is the exposure time, and $`\overline{r}`$ is the overall average count rate. For this work, we chose $`N=19`$ and the folding epoch to be MJD 52466.0 (which is roughly at the middle of the monitoring campaign). In the absence of any periodic or secular variations, $`S`$ follows the $`\chi ^2`$ distribution with $`N`$-$`1`$ degrees of freedom.

The true pulse period should correspond to where $`S`$ peaks. To accurately locate the peak, we fitted the $`SP`$ curve with a parabolic function in a small region around the peak. It should be noted that we did not subtract background from the light curves, in order to preserve simple statistics (Leahy et al. 1983). However, any non-Gaussian variation in the background rate would invalidate the basis for $`S`$ being $`\chi ^2`$ distributed. Fortunately, the effects of such variation seem to be quite small in our cases. A more serious complication arises from the fact that GX 1+4 itself exhibits aperiodic variations, as is evident by the presence of red noise in the power density spectrum (PDS; see Fig. 1 for an example). Therefore, $`S`$ is not a $`\chi _{N1}^2`$ random variable here, even in the absence of a periodic signal.

To reveal the underlying distribution of $`S`$, we constructed $`100`$ simulated light curves for each observation, following a recipe in Timmer & König (1995). In order to derive analytic PDS models for the simulations, we excluded from a measured PDS a few bins around the harmonics of the pulsed signal (with the exact number of bins removed depending on the exposure time of the observation), and subtracted off the Poisson noise (see, e.g., Zhang et al. 1996 for a description of the procedure). We then fitted the PDS continuum with a simple power law, which, by eye, seems to describe the data adequately. Note that in all cases except for one the short duration of the observation does not allow the sufficient averaging of the PDS, over multiple segments or frequency bins, to quantify the low-frequency noise with meaningful Gaussian errors. Consequently, we chose to seek best fits by following the least-squared algorithm, as opposed to the usual $`\chi ^2`$ minimization. For the long observation ($`15`$ ks), we binned the raw PDS by a factor of 10 and performed a $`\chi ^2`$ fit. The best-fit model is quite similiar to that from the least-squared fit and the reduced $`\chi ^2`$ of the fit is about 1.5.

The best-fit model was used to simulate light curves for each observation. Each simulated light curve was then re-randomized bin by bin (by replacing the number of counts, $`n`$, in a time bin with a random number drawn from a Poisson distribution with a mean value of $`n`$) to simulate noise due to photon counting processes. The resulted simulated light curves were then folded at the measured pulse period to yield a value of $`S`$. Finally, we constructed a histogram of $`S`$ for that observation (in the presence of red noise and photon counting noise alone). Fig. 2 shows the 68% confidence regions of $`S`$ for all observations, along with the values measured from actual data. The results are fairly robust with respect to modest amount of variation in the parameters of the input model. Clearly, the presence of red noise can seriously affect our ability in estimating the significance of the pulsed signal.

## 4 Results

Fig. 3 shows the flux history of GX 1+4 in 2002. The source was quite bright at the beginning of the year and was seen to make a transition to a faint state in mid-June. The average flux is roughly an order of magnitude lower in the faint state than in the bright state, although it varies significantly in both states. The minimum flux measured is about another order of magnitude lower, reaching $`3\times 10^{11}ergs\text{ }cm^2\text{ }s^1`$.

Toward the end of the year, GX 1+4 started to climb out of the faint state, as its timing and spectral properties gradually revert to those of the bright state (see discussions below). Unfortunately, our monitoring of the source was terminated too soon to cover the entire transition, because of the proximity of the Sun to the source.

### 4.1 Pulse Period and Pulse Fraction

X-ray pulsation is detected with high significance in all but four cases in the faint state. Fig. 4 (top panel) shows the $`S`$$`P`$ curve for one of the latter cases, which is basically flat. For comparison, we also show in the figure results for an adjacent observation (taken about a week earlier), where the pulsation is clearly detected. The measured pulse periods (from the PCA data) are also shown in Fig. 3. Note that no period is shown for the first observation, because it could not be reliably determined due to the very short exposure time of the observation (less than 4 pulse cycles). Although error bars are quite large in many cases (especially during the faint state), the spin-down trend can be seen from the figure. The best-fit spin-down rate is $`4.4\times 10^8\text{ }s/s`$, which is less than half of the BATSE average (Bildsten et al. 1997). The source reached a pulse period of $``$138.7 s, the largest that has ever been measured.

There are two sources that contribute to the uncertainty in measuring the pulse period: one is associated with the epoch-folding technique and the other is due to the presence of red noise. The former is well understood (see Larsson 1996), and we found that it was much smaller than the latter in all cases. We quantified the latter with simulations in the following manner. We added the observed pulsed signal to the simulated light curves (see § 3.2), for each of the observations where the pulsation was detected with high significance. We then derived pulse periods from the new simulated light curves through the same epoch-folding procedure. A histogram of pulse periods was constructed in the end and the 68% confidence region was determined. Therefore, the error bars have included the effects of intrinsic red noise of the source as well as of statistical fluctuations. We caution, however, that the reliability of the simulations depends sensitively on the accuracy of the input PDS model which we derived from the data.

The observed pulse profile varies significantly from observation to observation (i.e., over a timescale of days), much more so in the faint state. A thumbnail collection of pulse profiles are shown in Fig. 5. The profiles were derived by folding light curves (2–60 keV) with reference to the same epoch (MJD 52466.0), using the best-fit periods and period derivative (see Fig. 3). A feature that is common to a majority of the cases is the presence of a sharp dip in the pulse profile. The dips seem to occur at difference pulse phases. To be more certain about it, we experimented with deriving a phase coherence solution that connects the occurrence times of the dips, assuming they occur at a constant phase. In general, the pulse phase can be expressed as

$$\varphi (t)=\varphi (t_0)+\nu _0(tt_0)+\frac{1}{2}\dot{\nu _0}(tt_0)^2+\frac{1}{6}\ddot{\nu _0}(tt_0)^3+\mathrm{}$$

where $`t_0`$ is the reference epoch, and $`\nu _0`$, $`\dot{\nu _0}`$, and $`\ddot{\nu _0}`$ are the pulse frequency ($`1/P_0`$), its first and second derivatives at $`t=t_0`$, respectively. We started by breaking the longest observation (at MJD 52492) into six shorter segments. We folded the light curve of each segment and computed the time of the dip from the folded pulse profile. We did the same for two neighboring observations that were made about one week before and after the long observation. We were able to derive a timing solution (up to $`\ddot{\nu _0}`$) based on these data points. We then extended the time interval by adding one data point at a time and repeated the procedure. But, we failed to find a global phase coherent solution. Fig. 6 shows the residuals after the best-fit model has been subtracted from the data. This result seems to suggest that the dip does not always occur at a fixed pulse phase, although we cannot rule out the possibility that it is caused by non-vanishing higher-order period derivatives. Similar dips were seen previously (e.g., Dotani et al. 1989; Giles et al. 2000), and were interpreted as being either due to obstruction of the “hot spot” by the accretion column (e.g., Giles et al. 2000; Galloway et al. 2001) or to resonance scattering of photons in the accretion column above the magnetic pole (Dotani et al. 1989). However, it would seem difficult for either scenario to explain why the dip does not occur at a constant pulse phase.

As is obvious from Fig. 5, when GX 1+4 is relatively bright (the first 13 observations as well as the last few; see Fig. 3), its X-ray flux is strongly modulated. When the source reaches the faint state, the pulsed fraction is significantly reduced and the pulse profile takes on drastically different forms. In four of the cases (see Fig. 5), the folded light curves show no sign of periodic modulation. It should be kept in mind that GX 1+4 is known to exhibit red noise in its PDS (see Fig. 1 for an example). Consequently, much of the variation seen in the folded light curves can be caused by variations of the source that are associated with the red noise). To illustrate this point, we show, in Fig. 7, a comparison between a measured light curve and a simulated light curve for one of the observations in which no pulsation was detected. The non-Gaussian variation in the former can, therefore, be fully accounted for by the presence of red noise alone. We also investigated, with goodXenon data, the possibility that the signal might be detected only at higher energies for these cases. We constructed light curves that only include photons above roughly 7 keV and repeated the period-search procedure. No pulsation was detected. When the pulsation is not detected, the flux of GX 1+4 seems always to be $`2\times 10^{10}ergs\text{ }cm^2\text{ }s^1`$. On the other hand, in a few other cases, the flux of the source falls below that threshold yet the pulsation is detected with high significance. Therefore, the phenomenology is not simple.

To measure the amplitude of the pulsed signal, we defined the fractional peak-to-peak pulsed fraction as $`f_{pp}(F_{max}F_{min})/F_{max}`$, and computed it for each observation. The results are shown in Fig. 8. Of course, the pulsed fraction thus defined is not expected to be zero, even in the absence of any pulsed signal, because of statistical fluctuations and aperiodic variabilities. The results from simulations were again used to quantify the effects. Such “baseline” values of the pulsed fraction are shown in Fig. 8. Note that the pulsed fraction does not stand above the baseline in some cases, even though the pulsation is clearly detected. This is due to the combination of uncertainties in the simulation, small pulsed fraction, and poor statistics of the data.

### 4.2 Spectral Characteristics

The X-ray spectrum of GX 1+4 can be described by a model that consists of a cut-off power law and a Gaussian function. Also included in the model is the absorption along the line of sight. The observed continuum is, therefore, typical of accreting X-ray pulsars (White et al. 1983). The Gaussian component is needed to mimic the strong iron $`K_\alpha `$ line at $``$6.4 keV, which is known to be present in the X-ray spectrum of GX 1+4 (Kotani et al. 1999). This is the same model as the one that was adopted to study GX 1+4 during its 1996 transition to the faint state (Cui 1996) and thus facilitates direct comparisons of the results.

To account for residual calibration uncertainties, we followed the usual practice of including 1% systematic error in the analyses. The fits are quite satisfactory, with the reduced $`\chi ^2`$ values around $`1`$, except for a few cases at the beginning of the monitoring campaign when the source was bright. Even in the latter cases, the reduced $`\chi ^2`$ values are never more than $`1.7`$, with deviations residing mostly at the lowest energies, where the PCA calibration is less certain.

When it is bright, GX 1+4 shows little intrinsic spectral variability. The measured photon index falls in a narrow range of $`1.1`$$`1.3`$; the exponential cut-off energy remains roughly at $`8`$ keV; and the e-folding varies in the range of $`28`$$`40`$ keV. It is interesting to note that similar spectral characteristics were observed of the source during the last 3 observations, indicating the emergence of GX 1+4 from the faint state.

In contrast, the spectral properties of GX 1+4 appear to vary significantly in the faint state, on a timescale of days. Here, the photon index varies over a wide range of $`0.2`$$`2.0`$, with no apparent correlation with the measured flux. The spectrum of the source is the steepest in cases where the pulsation is unmeasurable. In most cases, the spectrum rolls over at energies between $`6`$$`13`$ keV, with the e-folding energy in the range of $`11`$$`45`$ keV. However, no spectral roll-over is observed when the pulsation ceases. To illustrate the effect, we made a composite spectrum by combining data from the four observations in which no pulsation is detected, in order to achieve sufficient statistics at high energies. Fig. 9 shows the spectrum, along with the best-fit model and the residuals. The presence of an iron $`K_\alpha `$ line is apparent; the continuum extends beyond 100 keV, following a simple power law with a photon index of roughly 1.6. To avoid divergence in overall luminosity, therefore, the spectrum must roll over at some energy, perhaps much beyond the HEXTE passing band. Very similar spectral variabilities were observed during the transition to a faint state in 1996 (Cui 1996).

The infered hydrogen column density varies in the ranges of $`4`$$`9`$ and $`3`$$`19\times 10^{22}\text{ }cm^2`$ for the bright state and faint state, respectively. It is not as well (sometimes, poorly) constrained in the latter case due to poorer statistics. These results are, again, comparable to those obtained previously (Cui 1996). The variability might be related to the changing physical conditions in the wind of the M giant companion star.

## 5 Discussion

In this work, we have presented evidence for the cessation of X-ray pulsation of GX 1+4, when the source reached a low-flux state. While it is not technically feasible to rule out the presence of weak pulsation in such a case, we can, at least, conclude that the X-ray pulsed fraction must have decreased dramatically as the source made a transition from the bright state to the faint state. Similar phenomenon was observed previously (Cui 1996), but we have carried out a more detailed analysis in work, taking into account the presence of intrinsic red noise with Monte-Carlo simulations. We expect that the derived error bars on the pulsed fraction measurements, as well as significance estimates regarding the detection of pulsation, are more realistic. This is a point worth emphasizing, because much of the variability seen in a folded light curve could be attributed to the aperiodic secular variations of the source.

One thing that we do notice in this work is that there does not appear to be a fixed flux threshold below which the pulsation ceases, as one might naively expect from the propeller effect. This complicates the interpretation. Perhaps, in the faint state, GX 1+4 was close to the transition threshold and it flipped and flopped stochastically between the propeller state and the non-propeller state. It might even be possible that there is some hysteresis associated with the transitions (see Fig. 2 of Lovelace, Romanova, & Bisnovatyi-Kogan 1999). If the propeller phenomenon is indeed observed, it allows a direct estimation of the magnetic field in GX 1+4.

### 5.1 Propeller Effect

For accreting pulsars, the determination of magnetospheric radius is still uncertain theoretically (e.g., Ghosh & Lamb 1979; Wang 1996). The geometry of the magnetic field is often assumed to consist of closed field lines, threading the accretion disk (Ghosh & Lamb 1979). More recently, other field configurations (with open field lines originating in the accretion disk) have been investigated and may, in fact, be important for the propeller effect (e.g., Lovelace, Romanova, & Bisnovatyi-Kogan 1999).

At the Alfvén radius, the magnetic pressure is balanced by the ram pressure of accreted matter. For a dipole field, the Alfvén radius is given by (Lamb, Pethick, & Pine 1973)

$$r_A=5.2\times 10^8\text{ }cm\text{ }L_{x,36}^{2/7}B_{12}^{4/7}M_{1.4}^{1/7}R_6^{10/7},$$
(2)
where $`L_{x,36}`$ is the bolometric X-ray luminosity in units of $`10^{36}\text{ }ergs\text{ }s^1`$, $`B_{12}`$ is the field strength at the surface of the neutron star in units of $`10^{12}\text{ }G`$, and $`M_{1.4}`$ and $`R_6`$ are the mass and radius of the neutron star in units of $`1.4M_{}`$ and $`10^6\text{ }cm`$, respectively. The magnetospheric radius is probably a fraction of the Alfvén radius, $`r_m=\eta r_A`$, where $`0<\eta 1`$ (e.g., Ghosh & Lamb 1979; Wang 1996; Lovelace, Romanova, & Bisnovatyi-Kogan 1999).

The co-rotation radius, $`r_{co}`$, is defined as where the Keplerian velocity is equal to the co-rotating velocity, $`(GM/r_{co})^{1/2}=\mathrm{\Omega }r_{co}`$, where $`\mathrm{\Omega }`$ is the angular velocity of the neutron star. Therefore, we have

$$r_{co}=4.5\times 10^9\text{ }cm\text{ }P_{138}^{2/3}M_{1.4}^{1/3},$$
(3)
where $`P_{138}`$ is the rotation period of the neutron star, in units of $`138\text{ }s`$. The propeller effect takes place when $`r_{co}=r_m`$. From equations (2) and (3), we find that the magnetic field strength is given by

$$B=4.4\times 10^{13}\text{ }G\text{ }\eta ^{7/4}P_{138}^{7/6}L_{x,36}^{1/2}M_{1.4}^{1/3}R_6^{5/2}.$$
(4)
For GX 1+4, the observed 2-60 keV X-ray fluxes are $`1.4`$$`1.7\times 10^{10}\text{ }erg\text{ }cm^2\text{ }s^1`$, for cases where the X-ray pulsation is undetectable. Assuming a distance of 6 kpc (which is probably accurate to within a factor of 2; Chakrabarty & Roche 1997), the observed luminosities are $`7\times 10^{35}\text{ }erg\text{ }s^1`$. Using this value in Eq. 4, we already have $`B>3.7\times 10^{13}\text{ }G`$. In some models, $`\eta `$ can be much less than 1 (e.g., $`0.3`$; Lovelace, Romanova, & Bisnovatyi-Kogan 1999), so the inferred magnetic field might be an order of magnitude larger. We must also make bolometric corrections, taking into account absorption at low energies and spectral roll-over at high energies, the latter of which is necessary to keep the integrated flux from diverging. While the low-energy correction is expected to be small ($`30\%`$), the high-energy correction can be quite large, depending on where the spectrum rolls over, given the relatively flat spectrum of the source. In extreme scenarios, therefore, the magnetic field of GX 1+4 could reach the magnetar regime!

### 5.2 Origin of the Persistent Emission in the Propeller State

In the propeller regime, the accreted matter can no longer be channeled to the surface of the neutron star by the magnetic field. Instead, it might be accumulated at the magnetosphere and gradually build up a dense envelope around it. In this scenario, the matter might be accumulated preferentially near the direction of the rotation axis and be expelled along the perpendicular direction (Wang & Robertson 1985). In the process, the matter could be heated up to tens of millions of degrees and produce thermal bremsstrahlung radiation. However, the spectrum of such radiation should roll over exponentially at tens of keV, which we did not observe in the case of GX 1+4. Therefore, the observed persistent emission is unlikely of magnetospheric origin.

The propeller effect was also invoked to explain the large luminosity variation of 4U 0115+63 near its pariastron passage (Campana et al. 2001). It was argued that the emission from the accretion disk would be important in the propeller regime. However, in the case of GX 1+4, the inner edge of the disk would be quite far away from the neutron star in the propeller regime and thus be too cool to contribute much to the observed X-ray emission. Moreover, the observed X-ray spectrum of GX 1+4 is clearly not that of disk emission.

On the other hand, we know that GX 1+4 has an M-giant companion star, which is expected to produce a relatively dense, slow stellar wind (Chakrabarty & Roche 1997). As the propelled matter plows into this wind at a velocity of the order of the escape velocity (Illarionov & Sunyaev 1975), $`(2GM/r_{co})^{1/2}2000`$ km/s, a strong shock might be formed. The shock could then produce a population of relativistic electrons via the first-order Fermi process. In this scenario, the observed persistent emission could simply be attributed to the synchrotron radiation of these electrons in the relatively strong magnetic field. The non-thermal nature of such emission is consistent with the lack of spectral roll-over in the propeller state. This emission mechanism is thought to be responsible for the unpulsed X–ray emission observed from the Be binary pulsar system PSR B1259-63 near periastron (Grove et al. 1995).

One potential complication might arise from the fact that GX 1+4 is located in the crowded Galactic center region. Since the flux of the source is quite low in the faint state, the issue of source confusion needs to be addressed. This can directly affect the derivation of pulsed fraction. More seriously, however, we need to be certain that GX 1+4 was detected at all in the faint state. We searched the existing X-ray catalogs for known sources within a $`1^{}`$ radius circle around GX 1+4, and only found a number of very faint ROSAT sources, some of which might be transient in nature. For a Crab-like spectrum, the fluxes of these sources are all below $`10^{13}\text{ }ergs\text{ }cm^2\text{ }s^1`$. This is consistent with results from our deep Chandra observation (with an exposure time of 60 ks) in the faint state: over the entire field of view ($`8.3`$′ $`\times `$ 50.6′ with the spectroscopic array), we detected only two very faint sources whose fluxes are more than 2 orders of magnitude lower than the faint-state flux of GX 1+4. Moreover, the fact that we did also detect pulsation at comparable fluxes is re-assuring. The presence of the characteristic iron $`K_\alpha `$ line in all observations adds another piece of evidence that the source never completely turned off its X-ray emission during our campaign.

We wish to thank Duncan Galloway for many helpful discussions and John Finley for his advice on pulse timing and comments on the manuscript. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. We gratefully acknowledge financial support from NASA through an LTSA grant (NAG5-9998) and a Chandra grant (GO2-3052X).


text:
# Dodecahedral space topology as an explanation for weak wide-angle temperature correlations in the cosmic microwave background

## 1 CMB Data Suggest Spherical Space with Dodecahedral Topology

Temperature fluctuations on the microwave sky may be expressed as a sum of spherical harmonics, just as music and other sounds may be expressed as a sum of ordinary harmonics. A musical note is the sum of a fundamental, a second harmonic, a third harmonic, and so on. The relative strengths of the harmonics – the note’s spectrum – determines the tone quality, distinguishing, say, a sustained middle C played on a flute from the same note played on a clarinet. Analogously, the temperature map on the microwave sky is the sum of spherical harmonics. The relative strengths of the harmonics – the power spectrum – is a signature of the physics and geometry of the universe. Indeed the power spectrum is the primary tool researchers use to test their models’ predictions against observed reality.

The infinite universe model gets into trouble at the low end of the power spectrum (Figure 1). The lowest harmonic – the dipole, with wave number $`\mathrm{}=1`$ – is unobservable because the Doppler effect of the solar system’s motion through space creates a dipole a hundred times stronger, swamping out the underlying cosmological dipole. The first observable harmonic is the quadrupole, with wave number $`\mathrm{}=2`$. WMAP found a quadrupole only about $`1/7`$ as strong as would be expected in an infinite flat space. The probability that this could happen by mere chance has been estimated at about a fifth of one percent (Spergel et al., 2003). The octopole term, with wave number $`\mathrm{}=3`$, is also weak at $`72\%`$ of the expected value, but not nearly so dramatic or significant as the quadrupole. For large values of $`\mathrm{}`$, ranging up to $`\mathrm{}=900`$ and corresponding to small-scale temperature fluctuations, the spectrum tracks the infinite universe predictions exceedingly well.

Cosmologists thus face the challenge of finding a model that accounts for the weak quadrupole while maintaining the success of the infinite flat universe model on small scales (high $`\mathrm{}`$). The weak wide-angle temperature correlations discussed in the introductory paragraph correspond directly to the weak quadrupole.

Microwave background temperature fluctuations arise primarily (but not exclusively) from density fluctuations in the early universe, because photons travelling from denser regions do a little extra work against gravity and therefore arrive cooler, while photons from less dense regions do less work against gravity and arrive warmer. The density fluctuations across space split into a sum of 3-dimensional harmonics – in effect the vibrational overtones of space itself – just as temperature fluctuations on the sky split into a sum of 2-dimensional spherical harmonics and a musical note splits into a sum of 1-dimensional harmonics. The low quadrupole implies a cut-off on the wavelengths of the 3-dimensional harmonics. Such a cut-off presents an awkward problem in infinite flat space, because it defines a preferred length scale in an otherwise scale-invariant space. A more natural explanation invokes a finite universe, where the size of space itself imposes a cut-off on the wavelengths (Figure 2). Just as the vibrations of a bell cannot be larger than the bell itself, the density fluctuations in space cannot be larger than space itself. While most potential spatial topologies fail to fit the WMAP results, the Poincaré dodecahedral space fits them strikingly well.

The Poincaré dodecahedral space is a dodecahedral block of space with opposite faces abstractly glued together, so objects passing out of the dodecahedron across any face return from the opposite face. Light travels across the faces in the same way, so if we sit inside the dodecahedron and look outward across a face, our line of sight re-enters the dodecahedron from the opposite face. We have the illusion of looking into an adjacent copy of the dodecahedron. If we take the original dodecahedral block of space not as a Euclidean dodecahedron (with edge angles $`117^{}`$) but as a spherical dodecahedron (with edge angles exactly $`120^{}`$), then adjacent images of the dodecahedron fit together snugly to tile the hypersphere (Figure 3b), analogously to the way adjacent images of spherical pentagons (with perfect $`120^{}`$ angles) fit snugly to tile an ordinary sphere (Figure 3a). Thus the Poincaré space is a positively curved space, with a multiply connected topology whose volume is 120 times smaller than that of the simply connected hypersphere (for a given curvature radius).

The Poincaré dodecahedral space’s power spectrum depends strongly on the assumed mass-energy density parameter $`\mathrm{\Omega }_0`$ (Figure 4). The octopole term ($`\mathrm{}=3`$) matches WMAP’s octopole best when $`1.010<\mathrm{\Omega }_0<1.014`$. Encouragingly, in the subinterval $`1.012<\mathrm{\Omega }_0<1.014`$ the quadrupole ($`\mathrm{}=2`$) also matches the WMAP value. More encouragingly still, this subinterval agrees well with observations, falling comfortably within WMAP’s best fit range of $`\mathrm{\Omega }_0=1.02\pm 0.02`$ (Bennett et al., 2003).

The excellent agreement with WMAP’s results is all the more striking because the Poincaré dodecahedral space offers no free parameters in its construction. The Poincaré space is rigid, meaning that geometrical considerations require a completely regular dodecahedron. By contrast, a 3-torus, which is nominally made by gluing opposite faces of a cube but may be freely deformed to any parallelepiped, has six degrees of freedom in its geometrical construction. Furthermore, the Poincaré space is globally homogeneous, meaning that its geometry - and therefore its power spectrum - looks statistically the same to all observers within it. By contrast a typical finite space looks different to observers sitting at different locations.

Confirmation of a positively curved universe ($`\mathrm{\Omega }_0>1`$) would require revisions to current theories of inflation, but the jury is still out on how severe those changes would be. Some researchers argue that positive curvature would not disrupt the overall mechanism and effects of inflation, but only limit the factor by which space expands during the inflationary epoch to about a factor of ten (Uzan, Kirchner and Ellis, 2003). Others claim that such models require fine-tuning and are less natural than the infinite flat space model (Linde, 2003).

Having accounted for the weak observed quadrupole, the Poincaré dodecahedral space will face two more experimental tests in the next few years:

* The Cornish-Spergel-Starkman circles–in–the–sky method (Cornish, Spergel and Starkman, 1998) predicts temperature correlations along matching circles in small multiconnected spaces such as this one. When $`\mathrm{\Omega }_01.013`$ the horizon radius is about 0.38 in units of the curvature radius, while the dodecahedron’s inradius and outradius are 0.31 and 0.39, respectively, in the same units; as a result, the volume of the physical space is only $`83\%`$ the volume of the horizon sphere. In this case the horizon sphere self-intersects in six pairs of circles of angular radius about $`35^{}`$, making the dodecahedral space a good candidate for circle detection if technical problems (galactic foreground removal, integrated Sachs-Wolfe effect, Doppler effect of plasma motion) can be overcome. Indeed the Poincaré dodecahedral space makes circle searching easier than in the general case, because the six pairs of matching circles must a priori lie in a symmetrical pattern like the faces of a dodecahedron, thus allowing the searcher to slightly relax the noise tolerances without increasing the danger of a false positive.
* The Poincaré dodecahedral space predicts $`\mathrm{\Omega }_01.013>1`$. The upcoming Planck surveyor data (or possibly even the existing WMAP data in conjunction with other data sets) should determine $`\mathrm{\Omega }_0`$ to within $`1\%`$. Finding $`\mathrm{\Omega }_0<1.01`$ would refute the Poincaré space as a cosmological model, while $`\mathrm{\Omega }_0>1.01`$ would provide strong evidence in its favour.

## 2 Conclusion

Since antiquity humans have wondered whether our universe is finite or infinite. For most of the past two millennia Europeans held the Aristotelian view of the universe as a finite ball with a spherical boundary. The invention of the telescope in 1608 revealed the universe to be much larger than Aristotle had imagined. Thus even though Galileo and Kepler clung to Aristotle s model, their successors Bruno, Descartes and especially Newton embraced the idea of infinite space. Nevertheless, some scientists were as uncomfortable with an infinite universe as they were with Aristotle’s hypothetical boundary. In 1854 Georg Riemann cut the Gordian knot by proposing the hypersphere as a model of a finite universe with no troublesome boundary. By 1890 Felix Klein discovered the more general concept of a multiconnected space (recall Figure 2) and during the early years of the twentieth century Einstein and others preferred finite universe models. Nevertheless, by the 1930’s the vast size of the observable universe had become known and the pendulum swung back towards infinite models. Now, after more than two millennia of speculation, observational data might finally settle this ancient question once and for all.

J.R.W. thanks the MacArthur Foundation for its generous support.