text:
# Effect of the galactic halo modeling on the DAMA/NaI annual modulation result: an extended analysis of the data for WIMPs with a purely spin-independent coupling

## I Introduction

In the last few years various technical approaches have been exploited in direct searches for Weakly Interacting Massive Particles (WIMP) direct\_exp . In particular, the DAMA/NaI Collaboration has collected a very large body of statistics which allows one to look for the distinctive signature in Dark Matter (DM) direct detection represented by the annual modulation of the rate damalast ; syst ; noi\_belli , an effect which is due to the rotation of the Earth around the Sun freese . The analysis of the DAMA/NaI data after 4 years of running damalast ; noi\_belli , corresponding to a total exposure of 57986 kg day, has indeed led to the observation of an annual–modulation effect, which does not appear to be related to any possible source of systematics syst . This exciting result has been analyzed under different hypotheses on the properties of WIMP dark matter: purely spin independent coupling damalast , mixed spin–coherent interaction dama\_mixed , inelastic dark matter dama\_inelastic . In the case of purely spin–independent interactions, the annual modulation result has been shown to be compatible with a galactic halo composed, at least partially, by relic neutralinos in different classes of supersymmetric models noi ; altri .

In the present analysis we consider the case of a WIMP with coherent interactions dominant over the spin–dependent ones. The analysis of the counting rate of any direct detection experiment may be done in terms of the WIMP mass $`m_W`$ and of the quantity $`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$, where $`\sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ is the WIMP–nucleon cross section for scalar interaction and $`\xi `$ 111In Ref. damalast  the same symbol indicates a different quantity: $`\xi =\rho _W/(0.3\mathrm{GeVcm}^3)`$. is the fractional amount of local non–baryonic DM density which is ascribed to the WIMP responsible for the effect ($`\xi 1`$) damalast ; noi . Performing a maximum–likelihood analysis of the data, the DAMA/NaI Collaboration has derived a region in the plane $`m_W`$$`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ which is compatible at 3–$`\sigma `$ C.L. to the observed effect of annual modulation. The properties of this region are sensitive to astrophysical inputs noi\_belli ; depart1 , some of which deserve a reanalysis and a deeper insight.

One of the main ingredients for the calculation of the expected rates is the distribution function of WIMPs in their six–dimensional phase space: $`F(\stackrel{}{r},\stackrel{}{v})d^3rd^3v`$ (where the position vector $`\stackrel{}{r}(x,y,z)`$ and the velocity vector $`\stackrel{}{v}(v_x,v_y,v_z)`$ are defined in the rest frame of the Galaxy). Direct detection rates $`R_{det}`$ depend on the distribution function (DF) at the Earth position in the Galaxy:

$$f(\stackrel{}{v})F(\stackrel{}{R}_0,\stackrel{}{v}),$$
(1)
where $`\stackrel{}{R}_0(R_0,0,0)`$ is the location of the Earth at a distance $`R_08.5`$ kpc from the galactic center and along the galactic plane. It is therefore clear that an accurate calculation of the expected detection rates requires a knowledge of the phase–space distribution function $`F(\stackrel{}{r},\stackrel{}{v})d^3rd^3v`$.

From the observational side, the most relevant piece of information coming from astrophysics is related to the rotational velocity of objects bounded to the Galaxy:

$$v_{rot}^2(r)=\frac{GM_{tot}(r)}{r},$$
(2)
where $`G`$ is the Newton’s constant and $`M_{tot}(r)`$ denotes the total mass contained inside the radius $`r|\stackrel{}{r}|`$:

$$M_{tot}(r)=_{r^{}<r}d^3r^{}\rho _{tot}(\stackrel{}{r^{}}).$$
(3)
The rotational velocity therefore depends on the total matter density distribution in space:

$$\rho _{tot}(\stackrel{}{r})=\rho _{DM}(\stackrel{}{r})+\rho _{vis}(\stackrel{}{r}),$$
(4)
where the DM density distribution function is:

$$\rho _{DM}(\stackrel{}{r})d^3vF(\stackrel{}{r},\stackrel{}{v}),$$
(5)
and $`\rho _{vis}(\stackrel{}{r})`$ represents the contribution to the matter density due to components other than the DM, like the disk and the bulge. The local values for the rotational velocity and for the DM matter density are denoted by $`v_0v_{rot}(\stackrel{}{R}_0)`$ and $`\rho _0\rho _{DM}(\stackrel{}{R}_0)`$ and they represent two key parameters in the calculation of WIMP direct detection rates, as it will be discussed in the following.

In order to calculate the DF of Eq. (1) one must invert Eq. (5) taking into account observational data. This problem is affected by degeneracies that cannot be solved without adding some piece of information. This explains why the velocity distribution represents one of the main source of uncertainty in the calculation of direct detection signals.

The usual approach to this problem consists in assuming that the system has some symmetry and that the distribution $`F`$ depends on the phase space parameters only through some integrals of motion (energy, angular momentum): this last condition automatically implies stationarity and that the Jeans’ equations are verified binney . The velocity ellipsoid $`\sigma _{ij}<v_iv_j>`$ may than be calculated as a function of the derivatives of the potential (which are related to the rotational velocity) by making use of Euler’s equation binney . Physically, this corresponds to impose hydrostatic equilibrium between pressure and gravitational attraction.

The most common and widely used example of such a procedure, and by far the simplest, is the isothermal sphere model. It consists in a spherical infinite system with a flat rotational curve, which automatically implies $`\rho (r)r^2`$ and the potential $`\mathrm{\Psi }\mathrm{ln}(r^2)`$ . The DF may be easily worked out, and turns out to be a Maxwellian: $`f(v)\mathrm{exp}(3v^2/(2v_{rms}^2))`$, where $`v|\stackrel{}{v}|`$ and $`v_{rms}`$ denotes the root mean squared velocity of the WIMPs. The isothermal sphere describes a self–gravitating gas of collisionless particles in thermal equilibrium, representing the highest entropy rearrangement of WIMPs in their phase space. A strong argument in favour of this last property is the “violent relaxation” model of Lynden–Bell violent\_relax , which indicates that the violently changing gravitational field of the newly formed Galaxy may have lead the non–interacting WIMPs to thermal equilibrium. Hydrostatic equilibrium and the assumption that the velocity ellipsoid is isotropic allows to calculate $`v_{rms}`$ through the relation: $`v_{rms}^2=3/2v_{rot}^2(R_0)`$. Due to its simplicity, the isothermal sphere model has become the “standard” assumption in the evaluation of DM expected rates, and has been used extensively in the literature, including the analysis of the DAMA/NaI modulation data damalast . However many of its underlying assumptions (sphericity of the halo, absence of rotation, isotropy of the dispersion tensor, flatness of the rotational curve) are not strongly constrained by astrophysical observations. Moreover the isothermal sphere is strictly unphysical and can only represent the behaviour of the inner part of physical systems, since it has a total infinite mass and needs some cut–off at large radii.

In light of the latest experimental data on WIMP direct searches, the issue of possible departures from the isothermal sphere model has gained interest and prompted several discussions noi\_belli ; depart1 ; anis . In the present paper we intend to analyze this issue in a systematic way, by employing a comprehensive set of self–consistent galactic halo models. Each model introduces a different degree of deviation from the simple isothermal sphere. We consider modifications in the velocity distribution function which are originated from a change of the gravitational potential or a change of the DM density profile evans\_pl ; evans\_log ; nfw ; moore ; kravtsov ; jaffe ; triaxial . We classify the different models depending on the symmetry properties of the galactic halo: spherical potential/density profile with an isotropic velocity dispersion; spherical potential/density profile with a non–isotropic velocity dispersion; axisymmetric models; triaxial models. For the axisymmetric models we also consider the possibility of having a co–rotating or counter–rotating halo.

The plan of the paper is as follows. In section II we introduce the formalism for the calculation of direct detection rates and the annual modulation signal, and summarize the procedure used by the DAMA/NaI Collaboration to determine the annual modulation region in the plane $`m_W`$$`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ for a purely spin–independent interacting WIMP. In Section III we describe the halo models that we intend to discuss and introduce a naming scheme that will be used throughout the paper. Section IV is devoted to the discussion of the constraints on the dark halo of our Galaxy coming from available observational data. In Section V the annual modulation region is calculated in a systematic way for all the models previously introduced, and the results are discussed. Finally, Section VI is devoted to our conclusions.

## II Direct detection rates and annual modulation effect

The expected differential event rate of a WIMP direct search experiment is given, for a monoatomic detector, by the expression:

$$\frac{dR_{det}}{dE_R}=N_T\frac{\rho _0}{m_W}𝑑\stackrel{}{w}f(\stackrel{}{w})w\frac{d\sigma }{dE_R}(w,E_R)$$
(6)
where $`N_T`$ is the number of the target nuclei per unit of mass, $`m_W`$ is the WIMP mass, $`\stackrel{}{w}`$ and $`f(\stackrel{}{w})`$ denote the WIMP velocity and DF in the Earth frame ($`w=|\stackrel{}{w}|`$), $`d\sigma /dE_R`$ is the WIMP–nucleus differential cross section and $`E_R`$ is the nuclear recoil energy. Notice that the detection rate is directly proportional to the local DM density $`\rho _0`$. The generalization of Eq.(6) to a diatomic detector, like NaI, is straightforward.

The differential cross section is, in general, a sum of a coherent and a spin–dependent contributions. In this paper we consider only the case of a WIMP whose interactions are dominated by the coherent part. In this case the rate may be expressed in terms of the WIMP–nucleon scalar cross section, $`\sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$, as:

$$\frac{d\sigma }{dE_R}\left(\frac{d\sigma }{dE_R}\right)_{coherent}\frac{F_n^2(E_R)}{E_R^{max}}\left(\frac{1+m_W/m_p}{1+m_W/m_N}\right)^2A^2\sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})},$$
(7)
where $`m_p`$ and $`m_N`$ are the proton and nucleus mass, $`A`$ is the nuclear mass number, $`E_R^{max}`$ is the maximal recoil energy and $`F_n(E_R)`$ is the nuclear form factor for coherent interactions which may be parametrized with the usual Helm expression helm .

The relation between the velocities $`\stackrel{}{v}`$ and $`\stackrel{}{w}`$ is given by:

$`\stackrel{}{v}`$ $`=\stackrel{}{v}_{}+\stackrel{}{w},`$ (8)
$`\stackrel{}{v}_{}`$ $`=\stackrel{}{v}_{}+\stackrel{}{v}_{rot},`$ (9)
where $`\stackrel{}{v}_{}`$ and $`\stackrel{}{v}_{}`$ denote the velocities of the Earth and the Sun in the Galactic rest frame ($`|\stackrel{}{v}_{}|v_0+`$12 km/sec), and $`\stackrel{}{v}_{rot}`$ is the Earth’s orbital velocity around the Sun ($`|\stackrel{}{v}_{rot}|=`$30 km/sec). Projecting Eq.(9) in the galactic plane, one gets:

$$|\stackrel{}{v}_{}|=|\stackrel{}{v}_{}|+|\stackrel{}{v}_{rot}|\mathrm{cos}\gamma \mathrm{cos}[\omega (tt_0)]$$
(10)
where $`\gamma `$ is the inclination of the plane of rotation with respect to the galactic one, $`\omega =2\pi /T`$ with $`T`$=1 year, and $`t_02^{nd}`$ June corresponds to the day when the Earth’s velocity is at its maximum.

The change of reference frame of Eqs. (8,9,10) introduces through the DF $`f(\stackrel{}{w})`$ a time dependence in the expected rate $`R_{det}`$. In order to exploit this time–dependence to extract the modulated signal from the measured counting rates, we follow the maximal–likelihood procedure of Ref. damalast , to which we refer for a detailed discussion. This procedure allows to determine the region in the plane ($`m_W`$$`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$) which is compatible with the modulation signal. A lower bound on $`m_W`$ at the value $`m_W=30`$ GeV is applied, to conform to the analysis of Ref.damalast ; noi\_belli .

The data we analyze in the present paper refer to the full set of data released so far by the DAMA/NaI Collaboration damalast  (DAMA/NaI 0–4), including also the upper limit on $`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ obtained by the same Collaboration, as discussed in Ref. damalast . The same values of the quenching factors and of the cut efficiences as in Refs. damalast ; noi\_belli  are used. We stress that also the determination of upper limits is affected by the choice of the WIMPs DF. This means that also when confronting upper limits one has to specify the galactic halo models which has been considered in the calculation. The formalism introduced in the present paper may in fact be used also to quantify the uncertainty in the determination of upper limits from direct detection experiments.

## III Halo models

Let us turn now to the discussion of the galactic halo models and of the techniques used to calculate the velocity distribution function $`f(\stackrel{}{v})`$. The different models are classified according to the symmetry properties of their matter density profile (or gravitational potential) and of the velocity distribution function. We define four classes:

* Spherically symmetric matter density with isotropic velocity dispersion;
* Spherically symmetric matter density with non–isotropic velocity dispersion;
* Axisymmetric models;
* Triaxial models.

All the models which we describe in this Section, and that will be used in the rest of the paper, are summarized in Table 2, where we also introduce a naming scheme that will be of practical use in the discussion.

### III.1 Spherically symmetric matter density with isotropic velocity dispersion

The first class of models is represented by those with a spherically symmetric matter density $`\rho (\stackrel{}{r})=\rho (r)`$ and isotropic velocity distribution $`f(\stackrel{}{v})=f(v)`$. These two conditions imply that the phase–space DF depends on the space and velocity variables only through the energy, which is an integral of motion: $`F(\stackrel{}{r},\stackrel{}{v})=F(ϵ)`$, where $`ϵ=\mathrm{\Psi }(r)v^2/2`$ is the relative energy (per unit mass) of the WIMP and $`\mathrm{\Psi }`$ is the relative potential, related to the total density $`\rho _{tot}`$ through the Poisson’s equation binney :

$$^2\mathrm{\Psi }=4\pi G\rho _{tot}.$$
(11)
Notice that $`\rho _{tot}`$ refers to all the matter components of the Galaxy, like the disk, the bulge or the halo, as written in Eq.(4).

Once the total potential $`\mathrm{\Psi }(r)`$ is known, the WIMP DF $`F(ϵ)`$ may be worked out by inverting Eq.(5). A change of variables from $`r`$ to $`\mathrm{\Psi }`$ allows to cast Eq.(5) as binney :

$$4\pi _0^\mathrm{\Psi }F(ϵ)\sqrt{2(\mathrm{\Psi }ϵ)}𝑑ϵ=\rho _{DM}(\mathrm{\Psi }),$$
(12)
By performing a Laplace inversion on Eq.(12), one obtains the Eddington formula binney :

$$F(ϵ)=\frac{1}{\sqrt{8}\pi ^2}\frac{d}{dϵ}_0^ϵ\frac{d\rho _{DM}(\mathrm{\Psi })}{d\mathrm{\Psi }}\frac{d\mathrm{\Psi }}{\sqrt{ϵ\mathrm{\Psi }}}.$$
(13)
In Eqs.(12,13) the normalization of $`\mathrm{\Psi }`$, which is defined through Poisson’s equation up to an arbitrary constant, is fixed by requiring that $`\mathrm{\Psi }(\mathrm{})=0`$. The velocity distribution function $`f(v)`$ which enters the calculation of direct detection rates is then obtained as in Eq.(1).

Eq.(13) shows that the dark matter DF $`F(\stackrel{}{r},\stackrel{}{v})`$ depends not only on the halo DM density distribution $`\rho _{DM}`$, but also on the density distribution of all the other galactic components. This has to be the case, since the DM particles move inside the gravitational potential generated by all the matter that makes up the Galaxy. A complete and rigorous determination of the DM DF will therefore require to model not only the galactic halo, but also all the other (disk, bulge) components of the Galaxy (see, for instance Ref. dehnen ).

However, WIMP direct detection is directly sensitive only to local properties of the Galaxy, and in particular the relevant parameters in the calculation of the detection rate are the local rotational velocity $`v_0`$ and the local DM matter density $`\rho _0`$. In fact, $`v_0`$ is directly related to the average WIMP kinetic energy, which is relevant in the scattering process with the nuclei of the detector, and the change of reference frame of Eqs. (8,9), which is crucial in determining the amount of annual modulation in the detection rate. Instead, $`\rho _0`$ is a sort of normalization factor for the direct detection rate $`R_{det}`$. The dependence of $`R_{det}`$ essentially on local parameters implies that a detailed modeling of the inner ($`rR_0`$) part of the Galaxy, where the disk and bulge components are more relevant and in general dominant over the halo, is not crucial for our analysis. Moreover, the matter density of non–halo components at the local position in the Galaxy ($`r=R_0`$) are no longer dominant with respect to the halo matter density at the same position (see, for instance Ref. dehnen ). The bulge, in fact, can be described by using a spheroidal density distribution which gives a sizeable contribution inside the first kpc from the galactic center, and it is truncated at about $`r`$ 2 kpc. The disk has an exponential distribution which in most of the models dies away at about 4 kpc from the galactic center. We will therefore assume in the following that in the outer Galaxy the dominant contribution to the matter density is given by the halo:

$$\rho _{vis}\rho _{DM}\mathrm{for}r>R_0.$$
(14)
The only basic information which is required from the non–halo components is their contribution to the local rotational velocity $`v_0`$:

$$v_0^2=v_{rot}^2(R_0)=\frac{G}{R_0}[M_{vis}+M_{halo}],$$
(15)
where

$$M_{vis(halo)}_{r^{}<R_0}d^3r^{}\rho _{vis(halo)}(r^{}).$$
(16)
A maximal halo occurs when $`M_{vis}M_{halo}`$: in this case almost all the local rotational velocity is supported by the halo and the local DM density $`\rho _0`$ gets its maximal value $`\rho _0^{max}`$ compatible with the given $`v_0`$. The opposite situation occurs when $`M_{vis}`$ assumes its maximal value compatible with observations: in this case, the local rotational velocity gets the maximal contribution from the non–halo components and only a fraction of $`v_0`$ is supported by the halo. At the same time $`\rho _0`$ gets its minimal value $`\rho _0^{min}`$, for the same $`v_0`$. The constraints on these parameters are discussed in Section IV.

From the point of view of calculating the DM DF (Eq.(13)), the occurrence of a maximal or non–maximal halo modifies the gravitational potential $`\mathrm{\Psi }(r)`$ and therefore the velocity distribution function $`f(v)`$ is affected. Indicating with $`\mathrm{\Psi }_0(r)`$ the potential for the maximal halo, the condition of Eq.(14) allows the generalization to the non–maximal case as:

$$\mathrm{\Psi }(r)=\frac{\rho _0}{\rho _0^{max}}\mathrm{\Psi }_0(r)+\left(1\frac{\rho _0}{\rho _0^{max}}\right)\frac{R_0}{r}v_0^2.$$
(17)
The condition of Eq. (14) allows to work out the total potential for the case of non-maximal halos without explicitly modeling the visible parts of the Galaxy (bulk, disk): all the dependence of these components are contained in $`M_{vis}`$.

As a comment, we notice that the presence of a non negligible contribution to $`v_0`$ from the non–halo components alters also the velocity distribution function of the isothermal sphere. The standard Maxwellian form for the isothermal sphere is in fact correct only for a maximal halo.

Now that we have discussed the procedure to calculate the velocity distribution function once the matter density of the DM is given, we proceed to introduce the different models.

The first type of model is a direct generalization of the isothermal sphere by introducing a core radius $`R_c`$ (model A1). The density profile is:

$$\rho _{DM}(r)=\frac{v_0^2}{4\pi G}\frac{3R_c^2+r^2}{(R_c^2+r^2)^2},$$
(18)
which corresponds to the following potential for a maximal halo:

$$\mathrm{\Psi }_0(r)=\frac{v_0^2}{2}\mathrm{ln}(R_c^2+r^2).$$
(19)
From the analytic form of this potential we will refer to this type of models as logarithmic models. The usual isothermal sphere (model A0) corresponds to the limit $`R_c0`$:

$$\rho _{DM}(r)=\frac{v_0^2}{4\pi G}\frac{1}{r^2},$$
(20)
which corresponds to the following potential for a maximal halo:

$$\mathrm{\Psi }_0(r)=\frac{v_0^2}{2}\mathrm{ln}(r^2).$$
(21)
For a maximal halo, the velocity distribution function arising from Eqs. (20,21) is the standard Maxwellian one. For a non–maximal halo, a deviation is present also in the $`R_c0`$, as discussed above. The rotational curve supported by $`\mathrm{\Psi }_0`$ of Eq. (19) is:

$$v_{rot}^2(r)=v_0^2\frac{r^2}{(R_c^2+r^2)}.$$
(22)
It is rising for small radii and at large radii becomes flat. In the $`R_C0`$ limit, $`v_{rot}=v_0`$ for all radii.

A second type of spherical models (A2 and A3) is defined by the following matter density evans\_pl :

$$\rho _{DM}(r)=\frac{\beta \mathrm{\Psi }_aR_c^\beta }{4\pi G}\frac{3R_c^2+r^2(1\beta )}{(R_c^2+r^2)^{(\beta +4)/2}}.$$
(23)
which corresponds to the following potential for a maximal halo:

$$\mathrm{\Psi }_0(r)=\frac{\mathrm{\Psi }_aR_c^\beta }{(R_c^2+r^2)^{\beta /2}}(\beta 0).$$
(24)
From the analytic form of this potential we will refer to this type of models as power–law models. They represent the spherical limit of the more general class of axisymmetric “power–law” models of Ref. evans\_pl  which will be discussed in Section III.2. The family of power–law models given by Eq.(24) is not defined for $`\beta =0`$. However, when $`\beta =0`$ is substituted in Eq. (23), the density of Eq. (18) is recovered. In fact the logarithmic model turns out to have the properties of the “missing” $`\beta =0`$ power law potential. Evaluating Eq.(23) for $`r=R_0`$ the parameter $`\mathrm{\Psi }_a`$ can be expressed in terms of the density $`\rho _0`$. The rotational velocity for the power–law models is given by:

$$v_{rot}^2=\frac{\beta \mathrm{\Psi }_aR_c^\beta r^2}{(R_c^2+r^2)^{(\beta +2)/2}},$$
(25)
and it is asymptotically falling with $`r`$ if $`\beta >0`$ (model A2) and rising if $`\beta <0`$ (model A3).

The last family of spherical models we consider (models A4–A7) is defined by the following matter density:

$$\rho _{DM}=\rho _0\left(\frac{R_0}{r}\right)^\gamma \left[\frac{1+(R_0/a)^\alpha }{1+(r/a)^\alpha }\right]^{(\beta \gamma )/\alpha }.$$
(26)
for the choice of values of the parameters $`\alpha `$, $`\beta `$, $`\gamma `$ and $`a`$ summarized in Table 1. Except for the Jaffe model (A4), the other three density profiles (A5, A6, A7) are obtained from numerical simulations of galaxy evolution.

### III.2 Spherically symmetric matter density with non–isotropic velocity dispersion

The procedure described in the previous Section can be generalized to the case of a non–isotropic velocity distribution, while keeping a spherically symmetric density profile. In this case, the most general DF is a function of $`ϵ`$ and of the magnitude of the angular momentum vector ($`L=|\stackrel{}{L}|`$): $`F=F(ϵ,L)`$. Among the different choices which can be performed, a particularly simple case is obtained when $`F`$ depends on $`ϵ`$ and $`L`$ only through the so called Osipkov–Merrit variable osipkov\_merrit ; binney :

$$Q=ϵ\frac{L^2}{2r_a^2},$$
(27)
where the parameter $`r_a`$ is related to the degree of anisotropy $`\beta _0`$ of the velocity dispersion tensor (evaluated at the Earth’s position) in the following way osipkov\_merrit :

$$\beta _0=1\frac{\overline{v}_\varphi ^2}{\overline{v}_r^2}=\frac{R_0^2}{R_0^2+r_a^2}.$$
(28)
Here the velocity is expressed in spherical coordinates and $`\overline{v}_\varphi =\overline{v}_\theta \overline{v}_r`$ (with $`\overline{v}_i^2<v_i^2><v_i>^2`$, $`i=r,\theta ,\varphi `$).

The corresponding DF can be obtained by solving a modified version of Eddington’s inversion formula, which is obtained by making the following substitutions in Eq.(13) binney :

$`ϵ`$ $``$ $`Q,`$ (29)
$`\rho _{DM}(r)`$ $``$ $`\rho _{Q,DM}(r)\left(1+{\displaystyle \frac{r^2}{r_a^2}}\right)\rho _{DM}(r).`$ (30)
The models we consider are the same as discussed in Section III.1: the logarithmic model of Eq. (18) (model B1), the power–law models of Eq. (23) (models B2 and B3) and the models defined by Eq.(26) (models B4–B7). The velocity distribution functions, which are obtained by solving the Eddington equation with the Osipkov–Merrit term, are therefore anisotropic with a degree of anisotropy controlled by the parameter $`\beta _0`$ related to $`r_a`$ as in Eq.(28).

### III.3 Axisymmetric models

In the case of axial symmetry the DF depends in general (at least) on two integrals of motion, the relative binding energy $`ϵ`$ and the component $`L_z`$ of the angular momentum along the axis of symmetry. The DF may be written quite generally as the sum of an even ($`F_+`$) and an odd ($`F_{}`$) contribution with respect to $`L_z`$:

$$F(ϵ,L_z)=F_+(ϵ,L_z)+F_{}(ϵ,L_z),$$
(31)
where:

$$F_\pm =\frac{1}{2}\left[F(ϵ,L_z)\pm F(ϵ,L_z)\right].$$
(32)
When Eq.(12) is extended to the axisymmetric case, the density $`\rho _{DM}`$ turns out to depend only on the even part $`F_+`$ binney , so that, by inverting it, the DF may be determined up to an arbitrary odd part $`F_{}`$. The problem of the determination of $`F_+`$ for an axisymmetric matter density is both analitycally and numerically hard to perform and actually it requires a double Laplace inversion on Eq.(12). However, for particular families of axisymmetric potentials this problem has been solved analytically by Evans evans\_pl ; evans\_log . These families are the axisymmetric generalization of the first two classes already introduced in Section III.1: the first family has a logarithmic potential, the second one has a power–law potential. We stress that these analytic solution for $`F_+`$ are obtained under the assumption that the halo potential is dominant over the other components. Therefore these solutions correspond to a maximal halo. Even though they do not represent the most general situation, their simplicity makes them of practical use and convenient for studying the axisymmetric case.

The first family of axisymmetric potential we consider is the logarithmic potential evans\_log  (models C1 and C2):

$$\mathrm{\Psi }_0(R,z)=\frac{v_0^2}{2}\mathrm{ln}(R_c^2+R^2+\frac{z^2}{q^2}),$$
(33)
where $`R^2=x^2+y^2`$ is the radial coordinate along the galactic plane, $`R_c`$ is the core radius and $`q`$ the flatness parameter. The corresponding DM density is:

$$\rho _{DM}(R,z)=\frac{v_0^2}{4\pi Gq^2}\frac{(2q^2+1)R_c^2+R^2+(2q^2)z^2}{(R_c^2+R^2+z^2q^2)^2}.$$
(34)
Eqs.(33) and (34) are the axisymmetric generalization of Eqs. (19) and (18). The corresponding rotational curve is obtained from Eq.(22) by substituting the radial coordinate $`r`$ with the radial coordinate in the galactic plane $`R`$.

By expressing $`z`$ as a function of $`\mathrm{\Psi }`$ through Eq.(33), the density of Eq.(34) may be decomposed as evans\_pl :

$$\rho =\rho _0(\mathrm{\Psi })+R^2\rho _1(\mathrm{\Psi }),$$
(35)
which allows to determine $`F_+`$ in the form:

$$F_+(ϵ,L_z)=F_+^0(ϵ)+L_z^2F_+^1(ϵ),$$
(36)
leading to a particularly simple analytic solution for the DF evans\_log . We give it for completeness in Appendix A. The relation of Eq.(35) is no longer valid for a non–maximal halo, since in that case, the change of variable from $`z`$ to $`\mathrm{\Psi }`$ is not determined by the potential of Eq.(33) but by an axisymmetric analogous of Eq.(17). Therefore the analytic expression given by Eq.(36) can be used only for maximal halos, i.e. $`\rho _0`$=$`\rho _0^{max}`$.

The second family of distribution functions is a generalization of the axisymmetric logarithmic potential to the case of an asymptotically non–flat rotational curve, while preserving the property of Eq. (35). This is obtained for the axisymmetric power–law potential evans\_pl  (models C3 and C4):

$$\mathrm{\Psi }_0(R,z)=\frac{\mathrm{\Psi }_aR_c^\beta }{(R_c^2+R^2+z^2q^2)^{\beta /2}}(\beta 0).$$
(37)
The corresponding matter density is:

$$\rho _{DM}(R,z)=\frac{\beta \mathrm{\Psi }_aR_c^\beta }{4\pi Gq^2}\frac{(2q^2+1)R_c^2+(1\beta q^2)R^2+[2q^2(1+\beta )]z^2}{(R_c^2+R^2+z^2q^2)^{(\beta +4)/2}}.$$
(38)
Evaluating Eq.(38) for $`R=R_0`$, $`z=0`$ the parameter $`\mathrm{\Psi }_a`$ can be expressed in terms of the density $`\rho _0`$. Eqs.(37) and (38) are the axisymmetric generalization of Eqs. (24) and (23) and they posses the same properties already discussed in Section III.1. As for the case of the logarithmic potential, the rotational curve is obtained from Eq.(25) by substituting the radial coordinate $`r`$ in Eq.(25) with the radial coordinate in the galactic plane $`R`$, and it is asymptotically falling with $`R`$ if $`\beta >`$0 (model C3) and rising if $`\beta <`$0 (model C4).

In analogy with the logarithmic case, also for the power–law model an analytic solution for the DF can be worked out evans\_pl  with the form of Eq. (36). This applies again only for a maximal halo. The analytic formulas for $`F_+`$ can be found for completeness in Appendix A.

#### III.3.1 Co–rotation and counter–rotation of the halo

As mentioned before, the DF for an axisymmetric model is known up to an arbitrary odd component $`F_{}`$. The DF we summarized above for the Evan’s models all refer to the pure even component: they all have $`F_{}(ϵ,L_Z)=0`$ and posses no bulk rotation. The case $`F_{}(ϵ,L_Z)0`$ corresponds to the case of a rotating halo, where the number of particles moving clockwise around the axis of symmetry is different from that in the opposite sense.

A family of DF’s with bulk rotation can be studied by constructing an explicit example for $`F_{}`$. This can be done, starting from a generic $`F_+`$, by considering the linear combination lynden\_corot ; kamion\_corot ; noi\_corot :

$$F_{}(ϵ,L_z)=F_{right}(ϵ,L_z)F_{left}(ϵ,L_z),$$
(39)
where

$$F_{right}(ϵ,L_z)=\{\begin{array}{cc}F_+(ϵ,L_z),\hfill & v_\varphi >0,\hfill \\ 0,\hfill & v_\varphi <0,\hfill \end{array}$$
(40)
and

$$F_{left}(ϵ,L_z)=\{\begin{array}{cc}0,\hfill & v_\varphi >0,\hfill \\ F_+(ϵ,L_z),\hfill & v_\varphi <0,\hfill \end{array}$$
(41)
The distributions $`F_{right}`$ and $`F_{left}`$ describe the configurations with maximal $`|\overline{v}_\varphi |`$ with the same density profile as $`F_+`$ lynden\_corot . A DF with an intermediate value of $`\overline{v}_\varphi `$ can be obtained as a linear combination of $`F_+`$ and $`F_{}`$, or, equivalently, of $`F_{left}`$ and $`F_{right}`$:

$$F(ϵ,L_z)=\eta F_{right}(ϵ,L_z)+(1\eta )F_{left}(ϵ,L_z).$$
(42)
The parameter $`\eta `$ ranges from 1 (maximal co–rotation) to 0 (maximal counter–rotation) and is related to the dimensionless spin parameter $`\lambda `$ of the Galaxy by: $`\lambda =0.36|\eta 0.5|`$ kamion\_corot . In order to be consistent with the available extensive numerical work on galaxy formation, $`\lambda `$ should not exceed the value 0.05 spingalaxy , implying $`0.36<\eta <0.64`$. For all the Evans model discussed in this Section, we will also study the co– and counter– rotating situations, adopting the two values $`\eta =0.36`$ and $`\eta =0.64`$

### III.4 Triaxial models

The last class of models we wish to discuss is represented by the triaxial potential discussed in Ref. triaxial :

$$\mathrm{\Psi }_0(x,y,z)=\frac{1}{2}v_0^2\mathrm{ln}\left(x^2+\frac{y^2}{p^2}+\frac{z^2}{q^2}\right),$$
(43)
which, for a maximal halo, corresponds to the DM density:

$$\rho _{DM}(x,y,z)=\frac{v_0^2}{4\pi G}\frac{Ax^2+By^2/p^2+Cz^2/q^2}{\left(x^2+y^2/p^2+z^2/q^2\right)^2}$$
(44)
where $`A=p^2+q^21`$, $`B=1+q^2p^2`$ and $`C=1+p^2q^2`$. In Ref. triaxial  the velocity DF $`f(\stackrel{}{v})`$ of the system is approximated by a triaxial Gaussian with semiaxes equal to the velocity dispersions as obtained by the solutions of the Jeans equations:

$`\overline{v}_r^2`$ $`=`$ $`{\displaystyle \frac{v_0^2}{(2+\delta )(p^2+q^21)}}`$ (45)
$`\overline{v}_\varphi ^2`$ $`=`$ $`{\displaystyle \frac{v_0^2(2q^21)}{2(p^2+q^21)}}`$ (46)
$`\overline{v}_\theta ^2`$ $`=`$ $`{\displaystyle \frac{v_0^2(2p^21)}{2(p^2+q^21)}}`$ (47)
when the Earth’s position is on the major axis of the equipotential ellipsoid (models D1 and D2), and:

$`\overline{v}_r^2`$ $`=`$ $`{\displaystyle \frac{v_0^2p^4}{(2+\delta )(1+q^2p^2)}}`$ (48)
$`\overline{v}_\varphi ^2`$ $`=`$ $`{\displaystyle \frac{v_0^2(2q^2p^2)}{2(1+q^2p^2)}}`$ (49)
$`\overline{v}_\theta ^2`$ $`=`$ $`{\displaystyle \frac{v_0^2(2p^2)}{2(1+q^2p^2)}}`$ (50)
when the Earth’s position is on the intermediate axis (models D3 and D4). In Eqs.(45,48) the quantity $`\delta `$ is a free parameter that in the spherical limit ($`p=q=1`$) quantifies the degree of anisotropy of the velocity dispersion tensor:

$$\frac{\overline{v}_\varphi ^2}{\overline{v}_r^2}=\frac{2+\delta }{2}.$$
(51)
## IV Constraining the models

Once a given model is chosen for the velocity distribution function of the dark matter particles, the parameters of the model have to be fixed using observational data. Unfortunately, due to its “darkness”, all our knowledge of the halo is of indirect nature turner ; dehnen  and it includes requirements on the circular rotational speed (constraints on its flatness and its value at the solar circle and in the outer regions of the Galaxy) as well as observational constraints on the local surface density of the disk and on the dispersion velocity of the bulge. In general, one should construct a composite model of the Galaxy where the DM is coupled to other components like the disk and the bulge, and the parameters describing the various components are varied independently requiring that the observational constraints are satisfied.

We have already noticed that WIMP direct detection rates are particularly sensitive to the value of the rotational velocity $`v_0`$ and the local DM density $`\rho _0`$ (both evaluated at the solar circle). The procedure we follow in order to determine the allowed ranges for $`v_0`$ and $`\rho _0`$ in each galactic model is explained in the following. First of all, the experimental informations we use are: the allowed range for the local rotational velocity, the amount of flatness of the rotational curve of our Galaxy and the maximal amount of non-halo components in the Galaxy. The first information directly fixes the allowed interval for $`v_0`$, irrespective of the galactic halo model. The other constraints are used in order to determine the allowed ranges for $`\rho _0`$, for each halo model.

The allowed interval for $`v_0`$ is:

$$v_0=(220\pm 50)\mathrm{km}\mathrm{sec}^1(90\%\mathrm{C}.\mathrm{L}.),$$
(52)
which conservatively relies on purely dynamical observations kochanek . Proper motion measurements of nearby stars cepheids  lead to similar estimates for the central value of $`v_0`$, with a significantly smaller uncertainty. However they are based on the assumption of circular orbit of the observed objects. For definiteness, we will use in the following three representative values for $`v_0`$, which correspond to its central value and to the boundaries of its allowed 90% C.L. range of Eq.(52): $`v_0=170,220,270`$ km sec-1.

For the three representative values of $`v_0`$ we then determine the corresponding allowed ranges for $`\rho _0`$. For each halo model and for each value of $`v_0`$, we calculate, as a function of $`\rho _0`$, two quantities: (i) the total amount of mass $`M_{vis}`$ in components other than the halo (e.g. : disk, bulge) which is necessary in order to match the given value of local rotational velocity $`v_0`$; (ii) the value of the rotational curve at a distance of 100 kpc from the center of the Galaxy: $`v_{rot}^{100}v_{rot}`$ ($`R=100`$ kpc). These two quantities are somewhat constrained from observations, even though their constraints are often obtained by using some degree of galactic modeling. We conservatively quote the following ranges turner ; dehnen :

$`110^{10}M_{}`$ $`<M_{vis}<`$ $`610^{10}M_{}`$ (53)
$`0.8v_0`$ $`<v_{rot}^{100}<`$ $`1.2v_0.`$ (54)
where $`M_{}`$ denotes the solar mass. The first constraint limits the amount of non-halo components, while the second is a constraints of “essentially flatness” rotational curve: only galactic halo models which provide a rotational curve which does not deviate from a flat one more than 20% at 100 kpc are accepted. The constraint of Eq.(54) is compatible with the estimates of the galactic mass at large radii as obtained by the dynamics of satellites of the Galaxy dehnen .

The behaviour of $`M_{vis}`$ and $`v_{rot}^{100}`$ as a function of $`\rho _0`$ is shown if Fig. 1 (which refers to $`v_0=220`$ km s-1), Fig. 2 ($`v_0=170`$ km s-1) and Fig. 3 ($`v_0=270`$ km s-1). Each line refers to a different halo model. The upper panels show that $`M_{vis}`$ is a decreasing function of $`\rho _0`$, since increasing the amount of dark matter in the Galaxy implies that less matter in other galactic components is required to support the rotational curve. On the other hand, the value of the rotational velocity in the outer Galaxy is totally supported by the dark halo, and it is larger for more massive haloes.

When the constraints expressed in Eqs. (53) and (54) are simultaneously applied, an allowed interval for $`\rho _0`$ may be derived for each halo model.

The procedure outlined above may be used as a simple recipe for identifying the intervals for the local density parameter. However, some caution must be taken in the application of the bounds on $`M_{vis}`$ in Eq.(53). As already discussed in the previous Section, for the models of class A and class B we solve the Eddington equation in order to determine the velocity distribution function of dark matter particles. In this case we can take into account also the situation in which the local rotational velocity is only partially supported by the halo, and therefore we can apply the limits to $`M_{vis}`$ given by Eq.(53). Instead, the analytic models of class C and class D may be applied only to the extreme case of a fully maximal–halo, since the analytic formulae for $`f(\stackrel{}{v})`$ are derived under the assumption that only the halo matter density is present. For these models (class C and class D) the only case we can deal with is that of a fully maximal halo, which corresponds to $`M_{vis}=0`$. This represents a conservative upper limit for $`\rho _0`$. From Fig. 1–3 we can see that upper bounds to $`\rho _0`$ obtained from the $`M_{vis}=0`$ limit are only a few percent larger than what is obtained by imposing the lower limit of Eq.(53): $`M_{vis}=110^{10}M_{}`$. For consistency and simplicity, we will use $`M_{vis}0`$ as a lower limit also for models of class A and B.

The allowed intervals for $`\rho _0`$ that we obtain by imposing the bounds on $`M_{vis}`$ and $`v_{rot}^{100}`$ are listed in Table 3. As discussed above, both values $`\rho _0^{min}`$ and $`\rho _0^{max}`$ will be used in the next Sections to perform the modulation analysis of the DAMA/NaI experimental data, while only $`\rho _0^{max}`$ will be used for the analysis of the data for models of class C and D. We stress that the reason for this stands in the fact that for the analytic models of classes C and D, we know $`f(\stackrel{}{v})`$ only for a maximal–halo.

## V Results and discussion

In this Section we make use of the halo models described in Section III, with the choice of parameters shown in Table 2, to analyze the annual–modulation signal present in the DAMA/NaI data damalast  in terms of relic WIMPs with purely coherent interactions. The procedure is the one outlined in Section II. The results are presented as 3$`\sigma `$ annual–modulation regions shown in the plane $`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ versus $`m_W`$. All figures are divided in three panels which correspond to the cases: $`v_0=170`$, 220, 270 km sec-1.

A general feature of all the models is that by raising the parameter $`v_0`$ the modulation region moves from the upper–right to the lower–left of the $`m_W`$$`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ plane. This is easily understood since, for a given DM density profile, higher values of $`v_0`$ imply higher values of $`\rho _0`$ (through Eq. (2)) and of the velocity ellipsoid $`\sigma _{ij}<v_iv_j>`$ (through the Jeans equations); the experimental value of the signal and the measured WIMP–nucleus recoil energy are fixed by the data, therefore the modulation region moves downward because:

$$\frac{dR_{det}}{dE_R}\rho _0\times \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}$$
(55)
(see Eq. (6)) and moves to lower masses because:

$$E_Rm_W\times <v^2>.$$
(56)
On the other hand, when one compares the different density profiles that we have discussed in the previous Sections, it is worth noticing that, for a given value of $`v_0`$, a stronger singularity in the galactic center lowers the value of $`\rho _0^{max}`$ (in order to keep constant the mass integral) with the consequence that the lower part of the modulation region rises. So the smallest values of $`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ are reached by the models with a less singular density profile.

Figures 4–18 show the result of the analysis for the models with a spherically symmetric density profile (models A0–7, B1–7). Models A0–7 have an isotropic velocity dispersion, while in models B1–7 a degree of anisotropy in the velocity dispersion is introduced through the Osipkov–Merrit term of Eq.(27). The amount of anisotropy is controlled by the $`\beta _0`$ parameter of Eq. (28) and it has been fixed at the value: $`\beta _0=0.4`$. This value corresponds to a radial anisotropy. In order to account for the uncertainty in the DM local density $`\rho _0`$, for each case two regions are given, the higher (lower) corresponding to $`\rho _0=\rho _0^{min}`$ ($`\rho _0^{max}`$), where $`\rho _0^{min}`$ and $`\rho _0^{max}`$ have been obtained as discussed in Section IV and are reported in Table 3. The numerical values of the parameters summarized in the third column of Table 2 have been chosen in order to ensure the compatibility of the curves of Figs. 1, 2, 3 with the constraints of Eqs.(53,54), discussed in Section IV. The Jaffe models A4, B4, which goes into the isothermal sphere when $`R_c\mathrm{}`$, have been calculated for the smallest allowed value of the core radius $`R_c`$ in order to examine the case of maximal departure from the usual scenario. As a consequence of this, for this model $`\rho _0^{min}\rho _0^{max}`$ and in Figs. 8,15 the upper and lower modulation regions are superimposed.

The effect of radial anisotropy in the velocity dispersion tensor ($`\overline{v}_r>\overline{v}_\theta =\overline{v}_\varphi `$), which occurs for the models of class B, may be seen by comparing Figs. 511 with the corresponding Figs. 1218. As a general feature, a reduction of the modulation effect is expected, since the WIMPs phase space is depopulated along the direction of the Sun’s velocity. This is confirmed by the fact that in most cases the modulation regions move upwards and widens, although the size of the effect can be small. The effect of radial anisotropy on the WIMP mass is more involved. In particular, the modulation regions for models B1, B3, B5, B6 extend to heavier WIMP masses compared to the corresponding isotropic cases, while for models B2, B4, B7 the region moves to smaller WIMP masses.

As already pointed out, for the models belonging to classes C and D only the regions for $`\rho _0`$=$`\rho _0^{max}`$ are shown. As far as the axisymmetric models of class C are concerned, they are shown in Fig. 19, 22, 25, 28. In each model the flatness parameter $`q`$ has been chosen in order to have the maximal flatness compatible with observations and with the positivity of the DF. The main effect induced by flatness is through the increase in the local density $`\rho _0`$, as can be seen in Table 3. As a consequence, the modulation regions for these models reach values of $`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ significantly below the levels of the spherical cases. In Figs. 20, 23, 26, 29 the same models are shown with a co–rotation effect of the halo implemented through Eq.(42) with $`\eta =0.64`$, while in Figs. 21, 24, 27, 30 for the same models a counter–rotation effect with $`\eta =0.36`$ has been introduced. The main consequence of halo co–rotation is a decrease of the relative velocity between WIMPs and the Earth. The energy of WIMPs is therefore smaller and in order to produce the same recoil energy in the detector the WIMPs have to be heavier. This can be verified in the figures, where the modulation region of co–rotating models may reach very high WIMP masses, even higher than 200 GeV. By the same token, in counter–rotating models the modulation region is shifted toward lower masses.

The peculiar shape of the modulation region of Fig. 20 deserves some comments. The two disconnected closed contours, which arise at different $`m_W`$ values, are indicative of the superposition in the WIMP phase space of two components with well separated r.m.s. velocities. This is exemplified in Fig. 35 for the models B1 ($`q=1`$) and C1 ($`q=1/\sqrt{2}`$), where the contour plots of the corresponding DF’s are plotted in the $`v`$$`v_\varphi `$ plane (in the galactic rest frame) and in the $`w`$$`w_\varphi `$ plane (in the Earth’s rest frame). It is evident from the figure that, in flattened models, small $`|L_z|`$ orbits are depopulated compared to the spherical case, leading to two well separated populations with $`v_\varphi >0`$ and $`v_\varphi <0`$. These two components have the same temperature in the galactic rest frame, but develop different r.m.s. velocities when boosted in the Earth’s rest frame. As shown in Fig. 35, the relative weight of the two populations can be tilted towards slow or fast WIMPs in the two cases of a co–rotating or counter–rotating halo, respectively. This explains why the co–rotating model of Fig. 20 develops a second minimum at high WIMP masses, and this is not observed in the corresponding non–rotating or counter–rotating cases. We have numerically verified that the peculiar disconnected region at high WIMP masses reduces in size when the flatness parameter is increased (since in this way the velocity distribution function becomes more similar to the non–flattened one), while it shifts towards lower masses when the core radius $`R_c`$ is decreased.

We conclude the discussion of our results with the triaxial models shown in Figs. 31, 32, 33, 34 where, to be definite, the same choice of parameters of Ref.triaxial  is adopted. For these models a general solution for the DF is not available. Only the velocity ellipsoid of Eq. (4548) is known, and it is used to fix the second moments of a non–isotropic Maxwellian. This explains why the shape of the modulation regions is quite similar to the standard case. In models D1 and D2 the Earth is assumed to be located on the major axis of the density ellipsoid, while in models D3 and D4 it is placed on the intermediate axis. Since in the two cases the Sun’s position $`R_0`$ is the same, $`\rho _0`$ is higher for models D1, D2 than for models D3, D4 (see Table 3). As a consequence, the modulation regions of Fig. 31, 32 reach smaller values of $`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ compared to those of Figs. 33, 34. Models D1 and D3 (D2 and D4) have $`\delta =1.78`$ ($`\delta =16`$), which implies a radial (tangential) anisotropy of the velocity ellipsoid (see Eq. (51)). Solving the Jeans equation for the potential of Eq.(43) in the spherical limit $`q=p=1`$ (which corresponds to a non-cored isothermal sphere) leads to the relation binney :

$$v_\varphi ^2+v_\theta ^2=\frac{3}{2}v_0^2.$$
(57)
Note that this property is a consequence of the flatness of the rotational curve, and would not be true, for instance, including a core radius in the potential. Eq.(57) implies that, in this model, the tangential components of the velocity ellipsoid are fixed by $`v_0`$. As a consequence of this, radial anisotropy ($`v_r>v_\theta =v_\varphi `$) corresponds to faster WIMPs and tangential anisotropy to slower WIMPs. This shifts the modulation regions towards smaller values of $`m_W`$ in Figs. 31, 33 and higher values of $`m_W`$ in Figs. 32, 34. This effect is sizeable in the case of tangential anisotropy, where the modulation region may extend up to $`m_W270`$ GeV.

The results of this section are summarized in Fig. 36, where all the modulation regions previously discussed have been plotted jointly. A convolution of all the regions may be indicative of the uncertainties in the determination of the WIMP modulation signal due to the modeling of the WIMP DF. As a final result, we show such a convolution in Fig. 37, where a single curve in the plane $`m_W`$$`\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}`$ is plotted by collecting the information contained in the analyses of all the non–rotating models considered in this paper. The region is compared with the original annual modulation region obtained in Ref. damalast  for an isothermal sphere model of the galactic halo with rotational velocity $`v_0=220`$ km s-1 and local dark matter density $`\rho _0=0.3`$ GeV cm-3. From Fig.37 we see that the DAMA/NaI annual modulation result is compatible with WIMPs masses up to $`m_W270`$ GeV and WIMP–nucleon cross sections in the interval: $`10^{10}\mathrm{nbarn}<\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}<6\times 10^8\mathrm{nbarn}`$ when the uncertainties in the WIMP velocity DF is taken into account. Co–rotating models with maximal corotation can extend the mass range even further, up to $`m_W`$ 500–900 GeV, for cross section of the order $`\mathrm{few}\times 10^9\mathrm{nbarn}<\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}<2\times 10^8\mathrm{nbarn}`$, as it can be seen, for instance, in Fig. 36.

## VI Conclusions

In the present paper we have extended a previuos analyses of the DAMA modulation experiment for the case of a WIMP with a purely spin-independent coupling, by discussing in detail the implications on the results of the uncertainties on the dark matter galactic velocity distribution. We have studied a large number of viable models which deviate from the standard isothermal sphere in the matter density profile, in the presence of anisotropies of the velocity dispersion tensor and in effects of rotation of the galactic halo. The different models have been classified according to the symmetry properties of their matter density profile (or gravitational potential) and of the velocity distribution function. We have specifically considered: (a) spherically symmetric matter density with isotropic velocity dispersion; (b) spherically symmetric matter density with non–isotropic velocity dispersion; (c) axisymmetric models; (d) triaxial models.

The different models have then been used to re–analyze the DAMA/NaI 0–4 data collected by the DAMA/NaI Collaboration damalast ; in particular a total exposure of 57986 kg day, which corresponds to 4 annual cycles, has led to the observation of an annual modulation effect. The hypothesis of WIMP annual modulation, already favoured in the previous studies damalast ; noi\_belli  by using an isothermal sphere, is confirmed in all the investigated scenarios, and the effects of the different halo models on the determination of the allowed maximum–likelihood region in the WIMP mass and WIMP–nucleon cross–section have been derived. We can summarize that the DAMA/NaI annual modulation result is compatible with WIMPs masses up to $`m_W270`$ GeV and WIMP–nucleon cross sections in the interval: $`10^{10}\mathrm{nbarn}<\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}<6\times 10^8\mathrm{nbarn}`$, when the uncertainties in the WIMP velocity DF is taken into account. When also co–rotation of the galactic halo is considered, the mass range extends further to $`m_W`$ 500–900 GeV, for cross section of the order $`\mathrm{few}\times 10^9\mathrm{nbarn}<\xi \sigma _{\mathrm{scalar}}^{(\mathrm{nucleon})}<2\times 10^8\mathrm{nbarn}`$. These intervals quantify the extent of the annual modulation region for WIMPs with purely spin–independent couplings, as due to uncertainties in the phase space distribution function of galactic WIMPs.

## Acknowledgements

We thank Prof. R. Bernabei and Prof. A. Bottino for useful discussions and suggestions. We are indebted with the DAMA/NaI Collaboration for allowing us to use its data for this investigation.

## Appendix A Axisymmetric models

In this appendix we give, for completeness, the analytic formulas of the DF’s for the axisymmetric potentials of Eqs. (33,37) adapted from Refs. evans\_log ; evans\_pl . All expressions are written in the reference frame of the Galaxy.

### A.1 Logarithmic potential

The DF for the logarithmic potential of Eq. (33) can be written as:

$$F(ϵ,L_z^2)=(A+B)\mathrm{exp}\left(\frac{2v^2}{v_0^2}\right)+C\mathrm{exp}\left(\frac{v^2}{v_0^2}\right),$$
(58)
where:

$`A`$ $`=`$ $`F_0\mathrm{\hspace{0.33em}4}\pi \left({\displaystyle \frac{2}{\pi }}\right)^{5/2}\left({\displaystyle \frac{v_\varphi }{v_0}}\right)^2{\displaystyle \frac{R_0^4}{(R_0^2+R_c^2)^2}}{\displaystyle \frac{1q^2}{q^2}},`$ (59)
$`B`$ $`=`$ $`F_0\mathrm{\hspace{0.33em}4}\pi \left({\displaystyle \frac{2}{\pi ^5}}\right)^{1/2}{\displaystyle \frac{R_0^2R_c^2}{(R_0^2+R_c^2)^2}}{\displaystyle \frac{1}{q^2}},`$ (60)
$`C`$ $`=`$ $`F_0\pi ^{3/2}{\displaystyle \frac{R_0^2}{R_0^2+R_c^2}}{\displaystyle \frac{2q^21}{q^2}},`$ (61)
and $`F_0=0.47`$ GeV cm$`{}_{}{}^{3}/v_0^3`$.

### A.2 Power–law potential

The DF for the power–law potential of Eq. (37) for $`\beta >0`$ can be written as:

$$F(ϵ,L_z^2)=A\stackrel{~}{ϵ}^{4/\beta 3/2}+B\stackrel{~}{ϵ}^{4/\beta 1/2}+C\stackrel{~}{ϵ}^{2/\beta 1/2},$$
(62)
where:

$$\stackrel{~}{ϵ}\frac{ϵ}{\mathrm{\Psi }_a}=\frac{\mathrm{\Psi }\frac{1}{2}v^2}{\mathrm{\Psi }_a}=\zeta ^\beta \frac{1}{2}\left(\frac{v}{v_1}\right)^2,$$
(63)
with $`\zeta =R_c/\sqrt{R_c^2+R_0^2}`$, while the velocity $`v_1=\sqrt{|\mathrm{\Psi }_a|}`$ is fixed through Eq. (38), and with suitable normalizations, can be cast in the form:

$$v_1=220\mathrm{km}\mathrm{sec}^1\left(\frac{\rho _0}{0.47\mathrm{GeVcm}^3}\right)^{\frac{1}{2}}\frac{R_0}{8.5\mathrm{kpc}}\frac{q(R_c^2+R_0^2)^{(\beta +4)/4}}{\sqrt{|\beta |}R_0R_c^{\beta /2}\sqrt{R_c^2(1+2q^2)+R_0^2(1\beta q^2)}},$$
(64)
The quantities $`A`$, $`B`$ and $`C`$ may be written as:

$`A`$ $`=`$ $`F_1\left({\displaystyle \frac{v_\varphi }{220\mathrm{km}\mathrm{sec}^1}}\right)^2\left({\displaystyle \frac{R_0}{R_c}}\right)^4{\displaystyle \frac{\mathrm{\Gamma }(2+4/\beta )}{2^{3/2}\pi ^{3/2}\mathrm{\Gamma }(4/\beta 1/2)}}\beta (\beta +2)({\displaystyle \frac{1}{q^2}}1),`$ (65)
$`B`$ $`=`$ $`F_1\left({\displaystyle \frac{v_1}{220\mathrm{km}\mathrm{sec}^1}}\right)^2\left({\displaystyle \frac{R_0}{R_c}}\right)^2{\displaystyle \frac{\mathrm{\Gamma }(2+4/\beta )}{2^{3/2}\pi ^{3/2}\mathrm{\Gamma }(1/2+4/\beta )}}\beta (\beta +2){\displaystyle \frac{1}{q^2}},`$ (66)
$`C`$ $`=`$ $`F_1\left({\displaystyle \frac{v_1}{220\mathrm{km}\mathrm{sec}^1}}\right)^2\left({\displaystyle \frac{R_0}{R_c}}\right)^2{\displaystyle \frac{\mathrm{\Gamma }(2+2/\beta )}{2^{3/2}\pi ^{3/2}\mathrm{\Gamma }(1/2+2/\beta )}}\beta \left[2{\displaystyle \frac{1+\beta }{q^2}}\right],`$ (67)
with $`F_1=0.47`$ GeV cm$`{}_{}{}^{3}/v_1^3`$.

For $`\beta <0`$ Eq. (62) still applies, with the following modifications (now $`\mathrm{\Psi }_a<0`$):

$$\stackrel{~}{ϵ}\frac{ϵ}{\mathrm{\Psi }_a}=\frac{\mathrm{\Psi }+\frac{1}{2}v^2}{\mathrm{\Psi }_a}=\zeta ^\beta \frac{1}{2}\left(\frac{v}{v_1}\right)^2,$$
(68)
where:

$`A`$ $`=`$ $`F_1\left({\displaystyle \frac{v_\varphi }{220\mathrm{km}\mathrm{sec}^1}}\right)^2\left({\displaystyle \frac{R_0}{R_c}}\right)^4{\displaystyle \frac{\mathrm{\Gamma }(3/24/\beta )}{2^{3/2}\pi ^{3/2}\mathrm{\Gamma }(14/\beta )}}\beta (\beta +2)(1{\displaystyle \frac{1}{q^2}}),`$ (69)
$`B`$ $`=`$ $`F_1\left({\displaystyle \frac{v_1}{220\mathrm{km}\mathrm{sec}^1}}\right)^2\left({\displaystyle \frac{R_0}{R_c}}\right)^2{\displaystyle \frac{\mathrm{\Gamma }(1/24/\beta )}{2^{3/2}\pi ^{3/2}\mathrm{\Gamma }(14/\beta )}}\beta (\beta +2){\displaystyle \frac{1}{q^2}},`$ (70)
$`C`$ $`=`$ $`F_1\left({\displaystyle \frac{v_1}{220\mathrm{km}\mathrm{sec}^1}}\right)^2\left({\displaystyle \frac{R_0}{R_c}}\right)^2{\displaystyle \frac{\mathrm{\Gamma }(1/22/\beta )}{2^{3/2}\pi ^{3/2}\mathrm{\Gamma }(12/\beta )}}\beta \left[{\displaystyle \frac{1+\beta }{q^2}}2\right],`$ (71)
with again $`F_1=0.47`$ GeV cm$`{}_{}{}^{3}/v_1^3`$.


format: text

text:
# On the limits of spectral methods for frequency estimation

## 1 Introduction

Nearly everywhere in science we encounter the problem of determining the frequency of oscillation of some autonomously oscillating system from which a signal has been measured over a finite time period. Most commonly, spectral methods are used: From the squared modulus of the Fourier transform of the signal (the periodogram) its expectation value (the power spectrum) is estimated, and the frequency is determined by some characterization of the position of a peak in the power spectrum (see, e.g., Mo et al., 1988; Muzi & Ebert, 1993; Timmer et al., 1996; Méndez et al., 1998; Godano & Capuano, 1999; Korhonen et al., 2002; Slavic et al., 2002). For linear systems driven by white noise it is easily seen that the periodogram or, equivalently, the empirical autocorrelation function do, in fact, contain all relevant information (Brockwell & Davis, 1991; Priestley, 1981).

For nonlinear oscillators, this is generally not the case. As was shown recently (Rossberg, 2002), the power spectrum of a noisy or chaotic (i.e., not perfectly periodic) nonlinear oscillator which oscillates with a given average frequency , i.e. the number of cycles per unit time, can be completely arbitrary. In particular, the power spectrum is unrelated to the frequency. The number of oscillation cycles is here defined by the number of transition of the trajectory of the signal in delay space through a Poincaré section. Since this corresponds to a purely topological relation between the trajectory and the Poincaré section, the corresponding measure of frequency was named *topological frequency* $`\omega _{\text{top}}`$. Of course, intuition and experience suggest that, for typical experimental time series, the result of using some standard method to determine the oscillation frequency from the power spectrum would not differ too much from the topological frequency. But from the power spectrum alone no upper limit for this difference can be derived. Nor seems any other efficient method to be known which would allow to determine an upper limit for this difference.

After giving some examples for oscillations which exhibit this unrelatedness of topological frequency and power spectrum in Section 2, a result is presented which is in some sense stronger than the result cited above. A precise determination of frequency is not only impossible from the power spectrum of a time series, but also from the periodogram. Since the power spectrum is, in one way or another, estimated from the periodogram by an averaging procedure, the estimated power spectrum contains less information than the periodogram. One might conjecture that from the additional information contained in the periodogram, such as higher order correlations between the amplitudes of individual Fourier modes, frequency information can be extracted. As is shown in Section 3, this is not the case: Examples for time series are given which obviously have different topological frequencies, but, by construction, *identical* periodograms. The use of the periodograms of the raw data is justified for spectral estimation only if end effects can be neglected, for example, when the time series are long enough. In Section 4 it is shown that, in fact, for long enough time series the result holds valid also when end effects are explicitly taken account by spectral estimators.

## 2 Oscillation frequencies and power spectra of some nonlinear oscillators

As a first example, consider the noisy, weakly-nonlinear oscillator described by a complex amplitude $`A(t)`$ with dynamics given by the noisy Landau-Stuart equation (Risken, 1989)

$`\dot{A}=(ϵ+i\omega _0)A(1+ig_i)|A|^2A+\eta (t),`$ (1)
where $`ϵ`$, $`\omega _0`$, and $`g_i`$ are real and $`\eta (t)`$ denotes complex, white noise with correlations

$`\eta (t)\eta (t^{})=0,\eta (t)\eta (t^{})^{}=4\delta (tt^{})`$ (2)
\[ $``$ complex conjugation, $``$ $``$ expectation value\]. In a certain sense (Arnold, 1998), this system universally describes noisy oscillations in the vicinity of a Hopf bifurcation. In general (i.e., with $`g_i0`$) the linear frequency $`\omega _0`$, the *spectral peak frequency* $`\omega _{\text{peak}}`$, i.e. the position of the maximum of the power spectrum $`S_A(\omega )`$, the average frequency or *phase frequency*

$`\omega _{\text{ph},A}:=\omega _i,\omega _i:=\mathrm{Im}\left\{\dot{A}/A\right\},`$ (3)
and the *mean frequency* of $`A(t)`$

$`\omega _{\text{mean},A}:={\displaystyle \frac{\omega _i|A|^2}{|A|^2}}={\displaystyle \frac{\mathrm{Im}\dot{A}A^{}}{|A|^2}}={\displaystyle \frac{\omega S_A(\omega )𝑑\omega }{S_A(\omega )𝑑\omega }},`$ (4)
are all different; see Fig. 1. \[Definitions (3,4) are sometimes restricted to *analytic signals*, characterized by $`S_A(\omega )=0`$ for $`\omega 0`$, which are in one-to-one correspondence with the zero-mean, real-valued signals $`\mathrm{Re}\{A(t)\}`$. See Boashash (1992) for the history.\]

The phase frequency measures the average number of circulations around the point $`A=0`$ in phase space per unit time (decompose $`A(t)=a(t)e^{i\varphi (t)}`$ to see this). Thus, with the natural assumption that a circulation around the origin of phase space corresponds to a cycle of oscillation – which is justified by symmetry considerations – the phase frequency measures cycles per unit time, i.e., the oscillation frequency. Although this concept of an oscillation cycle is slightly different from that used for the definition of topological frequency, both frequency measures often give identical results, in particular when the oscillation amplitude does not come too close to zero. For a discussion of the effect of filtering on the phase frequency, see Rossberg (2002). In contrast to the phase frequency, the spectral peak frequency $`\omega _{\text{peak},A}`$ and the mean frequency $`\omega _{\text{mean},A}`$ are frequency measures which can be obtained directly from the power spectrum. As can be seen from Eqs. (3) and (4), $`\omega _{\text{ph},A}=\omega _{\text{mean},A}`$ if and only if $`\omega _i`$ and $`|A|^2`$ are uncorrelated. This is always the case for linear, noise-driven oscillators, but for general nonlinear oscillators it is not.

When it is known that the dynamics of $`A(t)`$ is of the form of Eqs. (1) and (2) and the power spectrum $`S_A(\omega )`$ of $`A(t)`$ is precisely determined, then it is of course possible to obtain the free parameters of Eq. (1) from a fitting procedure and to use these to determine $`\omega _{\text{ph},A}`$. But none of these assumptions is typically satisfied in practice. Equations (1) and (2) are exact only in a particular limit, and $`S_A(\omega )`$ is usually biased by some nontrivial transfer function inherent in the measurement process. Then the information about $`\omega _{\text{ph},A}`$ cannot be extracted from the power spectrum anymore. The situation becomes aggravated by the fact that, as was shown by Seybold & Risken (1974), $`S_A(\omega )`$ is always well approximated by a Lorenzian, uniformly for all $`ϵ`$.

As a second example, an experimental time series from the hand tremor of a subject diagnosed with Parkinson’s disease shall be examined. The data was recorded by an accelerometer attached to the hand. It consists of $`10240`$ samples recorded at a rate of $`300\mathrm{Hz}`$. The data is available on the internet222http://webber.physik.uni-freiburg.de/$``$jeti/tremordaten/park/0308-2-ali.dat.

Figure 2a shows part of the time series. In order to obtain an estimate of the power spectrum (Fig. 2b), the periodogram of the (tapered) time series was smoothed using a window of width $`h=10`$ with weights $`w(i)=\mathrm{max}(0,1i^2/h^2)`$. Its peak frequency is at $`\omega _{\text{peak}}/2\pi =9.49\mathrm{Hz}`$, corresponding to 324 oscillations over the time series. In order to determine the phase frequency, the (band-pass filtered) analytic signal corresponding to the time series is obtained by convoluting it with the Morlet wavelet

$`w(t)={\displaystyle \frac{1}{\pi }}\mathrm{exp}\left(i\omega _0t{\displaystyle \frac{1}{2}}(\mathrm{\Delta }\omega t)^2\right)`$ (5)
with $`\omega _0=\omega _{\text{peak}}`$ and $`\mathrm{\Delta }\omega =\omega _{\text{peak}}/6`$. The resulting phase frequency $`\omega _{ph}=9.58\mathrm{Hz}`$ (Fig. 2b) is robust under large variations of the filter width $`\mathrm{\Delta }\omega `$ and mid frequency $`\omega _0`$. It is also identical to the topological frequency obtained by band-pass filtering the time series by a convolution with $`\mathrm{Re}\{w(t)\}`$ and a successive 2D delay embedding, e.g., with delay $`\tau =(23/300)\mathrm{s}`$. Unfortunately, the investigated time series is too short to yield the observed difference between peak frequency and topological frequency significant in a statistical sense, as a careful statistical analysis using the method of Timmer et al. (1997) shows: The $`95\%`$ confidence interval for the peak frequency extends from $`9.40\mathrm{Hz}`$ to $`9.61\mathrm{Hz}`$. It would be desirable to repeat the analysis on longer samples.

It should be emphasized that the point of Rossberg (2002) is that the topological or phase frequency cannot be derived from the power spectrum *as a whole*. Determination of the spectral peak frequency is only the most obvious attempt to do so. The failure of this particular approach has been noticed before, see, e.g., Pikovsky et al. (1997).

## 3 Does the periodogram help?

Is it possible to identify a relative shift between spectral peak frequency and topological or phase frequency by using a different method to evaluate the periodogram? As mentioned before, the answer to be given here is negative: there is *no* general way to tell from the periodograms (or the empirical autocorrelation functions) of time series that they oscillate with different frequencies.

For simplicity, only real valued time series $`x(t)`$ shall be considered here for which the number of oscillations in a given time interval of length $`T`$ can be determined by a simple 2D delay embedding. If there is a delay $`\mathrm{\Delta }t`$, such that the trajectory of the delay vector $`(x(t),x(t\mathrm{\Delta }t))`$ encircles the origin of the coordinate system at finite distance, the number of circulations $`m`$ gives the number of oscillations in the time series. The resulting frequency is then $`2\pi m/T`$. This is the definition of topological frequency restricted to 2D embedding. It is very similar to the one that can be obtained from of the definition of instantaneous phase $`\varphi (t)`$ by Pikovsky et al. (1997) when defining $`\omega `$ as the temporal average of $`d\varphi /dt`$. However, their definition is based on the actual phase space of a dynamical system and not on its approximate reconstruction from a time series by delay embedding. This distinction becomes particularly relevant when random time series, as they are modeled by the method described here, are evaluated.

We shall now outline an iterative algorithm that randomly generates pairs of time series $`x_1(t)`$, $`x_2(t)`$ of the type described above. The two time series generated by this algorithm have identical periodograms but differ in the number of oscillation cycles. (The algorithm can also be considered as a source of random time series $`X_1(t)`$, $`X_2(t)`$ which have identical power spectra but different expectation values for the oscillation frequencies – in fact, the frequencies are fixed). The algorithm uses some ideas from an algorithm of Schreiber & Schmitz (1996) that generates time series with periodograms nearly identical and distributions of values identical to a given times series.

The main part of our algorithm generates two complex-valued random time series $`y_1(t)`$, $`y_2(t)`$ with $`N`$ samples at $`t=0,\mathrm{},N1`$. When continued periodically \[$`y_l(t)=y_l(t+N)`$\], the distribution of the output is invariant under translation in time. The periodograms of the two time series are identical. The moduli $`|y_l(t)|`$ follow a predetermined narrow distribution with values near one, i.e., the values of $`y_l(t)`$ lie near the unit circle in the complex plane. $`y_1(t)`$ and $`y_2(t)`$ encircle the origin of the complex plane $`k_1`$ and $`k_2`$ ($`k_1k_2`$) times respectively. Using these time series, the oscillatory, real-valued time series $`x_1(t)`$, $`x_2(t)`$ ($`x_l(t)=\mathrm{Re}\{\mathrm{exp}(i\omega _0t)y_l(t)`$) with identical periodograms are constructed.

In order to obtain the final time series $`y_l(t)`$, the following two basic steps are performed iterative, starting with two time series which encircle the origin of the complex plane $`k_1`$ respectively $`k_2`$ ($`k_1k_2`$) times, but do not yet have identical periodograms: In Fourier space, the magnitudes of the Fourier modes of the two time series are adjusted, such that the magnitudes of the Fourier modes at the same frequencies approach each other, while the phases of the Fourier modes are kept constant (Steps 4 and 5 in the detailed description below). Then, back in physical “space”, i.e., using the actual time series, a nonlinear, monotonously increasing transformation of the magnitudes of the values for each time is applied in such a way that the desired distribution of the complex time series near the unit circle (see above) is enforced (Step 10 below). If the phases were kept constant in this second basic step, the algorithm would (as numerical experiments show) typically converge to a trivial solution with $`y_1(t)=y_2(t)`$. Obviously, the number of circulations of the origin of the complex plane would then have changed with respect to the initial condition – at least for one of the two time series. In order to avoid these phase slips, the difference in phase angle between successive values is, for both time series, enforced to lie within a certain, small rage, by adjusting, at each iteration, also the phase angles (Steps 1,6-9). The details of these manipulations of the phases, as well as the precise choice of the initial time series (Steps 2 and 3), contain random elements, as will be explained below. Experience shows that convergence of the algorithm is, although not guaranteed, with appropriate choices of the tuning parameters of the algorithm (see below), highly probable.

This is a step-by-step description of this algorithm. We give the values of tuning parameters as they were used to generate the sample Pair A shown in Fig. (3), which consists of two time-series with $`N=256`$ points, in parenthesis.

1. Generate two random permutations $`\mathrm{\Pi }_l:(0\mathrm{}(N1))(0\mathrm{}(N1))`$ ($`l=1,2`$) of length $`N`$, with equal probability $`1/N!`$ for each permutation, to be used later.
2. Initialize the series $`y_1(t)`$ and $`y_2(t)`$ as $`y_l(t)=\mathrm{exp}(2\pi ik_lt/N)`$ at $`t=0,\mathrm{},N1`$ with small integer $`k_l`$ (e.g., $`k_1=0`$ and $`k_2=1`$).
3. In order to excite also the remaining Fourier modes, add some noise to $`y_1(t)`$ and $`y_2(t)`$ (e.g., complex Gaussian with variance $`2\times 10^4`$).

Now the main loop (steps 4-10) starts.
4. Calculate the discrete Fourier transforms $`\stackrel{~}{y}_l(\omega )`$ of the time series $`y_l(t)`$. Increase or decrease the values of $`|\stackrel{~}{y}_l(\omega )|`$ such that the distance between $`|\stackrel{~}{y}_1(\omega )|`$ and $`|\stackrel{~}{y}_2(\omega )|`$ decreases by a fixed percentage (e.g., $`20\%`$) while keeping $`\mathrm{arg}(\stackrel{~}{y}_l(\omega ))`$ and $`|\stackrel{~}{y}_1(\omega )|+|\stackrel{~}{y}_2(\omega )|`$ fixed. In this step the two periodograms are forced to approach each other. Notice, however, that the manipulations in steps 6-10 may have an adversary effect. In order to have the algorithm as a whole converge, some balancing of the tuning parameters is therefore required.
5. Obtain the inverse Fourier transforms of the results and assign them to $`y_1(t)`$ and $`y_2(t)`$.
6. Calculate the increments of $`\mathrm{arg}(y_l(t))`$ between successive samples. Apply the permutations $`\mathrm{\Pi }_l`$ to the sequences of phase increments (see below for an explanation).
7. For phase increments of absolute value larger than some threshold (e.g., 0.4) distribute part of the increment (e.g., $`1\%`$) symmetrically over the preceeding and succeeding increments such that the sum of the three remains unchanged. Do this step with parallel update for all increments.
8. Repeat step 7 until no phase increment exceeds the threshold.
9. Obtain the inverse permutations $`\mathrm{\Pi }_l^1`$ of these phase increments, and sum them up to obtain two new time series for the phases $`\mathrm{arg}(y_l(t))`$.
10. For each series, find a polynomial $`P_l()`$ of low order (e.g., 9th order), such that the distribution of $`P_l(|y_l(t)|)`$ fits best with the desired distribution of the arguments (we use a Gaussian with mean $`m=1`$ and variance $`\sigma ^2=0.04`$).333$`P_l()`$ is determined such that    $`\mathrm{min}={\displaystyle \underset{s=0}{\overset{N1}{}}}\left[P_l\left(\left|y_l\left(Q_l(s)\right)\right|\right)m2^{1/2}\sigma \stackrel{1}{\text{erf}}\left({\displaystyle \frac{2s+1}{N}}1\right)\right]^2,`$  (6)  where $`Q_l()`$ is the permutation which rank orders the sequence $`\{|y_l(t)|\}_t`$ and $`\text{erf}^1()`$ denotes the inverse error function given by    $`x={\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle _0^{\text{erf}^1(x)}}e^{u^2}𝑑u\text{for all }1<x<1\text{.}`$  (7)  <sup>,</sup><sup>4</sup><sup>4</sup>4Schreiber & Schmitz (1996) achieve a similar result by using, instead of $`P_l()`$, functions $`M()`$ which map the sets of values of the time series (e.g.,$`\{|y_1(t)|\}_t`$) onto a given target set of values in such a way that rank ordering is preserved. The method used here has the advantage over their method that the periodograms can be matched to arbitrary accuracy. Construct two new complex-valued time series by using the arguments $`P_l(|y_l(t)|)`$ and the phases from step 9, and assign them to $`y_1(t)`$ and $`y_2(t)`$.
11. Repeat steps 410 until the periodograms of $`y_1(t)`$ and $`y_1(t)`$ converge to the same values, i.e., $`\delta :=N^1_\omega (|\stackrel{~}{y}_1(\omega )||\stackrel{~}{y}_2(\omega )|)^2<ϵ`$, with some small $`ϵ>0`$ (e.g., $`ϵ=10^6`$). For larger $`N`$ convergence is substantially improved by judiciously choosing new permutations $`\mathrm{\Pi }_l`$ from time to time (see below).
12. Multiply random phases $`\mathrm{exp}(i\varphi _l)`$ to $`y_1(t)`$ and $`y_2(t)`$ in order to make the distributions of $`x_1(t)`$ and $`x_2(t)`$ obtained in the next step time-translation invariant.
13. Calculate the real-valued time series $`x_l(t)=\mathrm{Re}\{\mathrm{exp}(i\omega _0t)y_l(t)\}`$ with some not too small $`\omega _0`$ (e.g., $`\omega _0=5\pi /32`$).
14. As in steps 4 and 5, remove the small differences in the periodograms of $`x_1(t)`$ and $`x_2(t)`$ resulting from the $`\mathrm{Re}\{\}`$ operation in step 13. But now, do it in a single ($`100\%`$) step.

A typical result of running this algorithm is shown in Fig. 3. With a high enough frequency of the fast modulation $`\omega _0`$, the algorithm guarantees that the origin of delay space \[$`\mathrm{\Delta }t\pi /(2\omega _0)`$\] is encircled at finite distance (see Fig. 4).

Now, some details concerning the manipulation of phases (Steps 6-9) shall be explained. The permutations of the phase increments prior to the relaxation of large phase jumps in steps 7 and 8 prevent the algorithm from running into solutions with localized regions of fast drift of phase, which do not occur in typical experimental time series. However, it turns out that, with permutations selected with equal probability from the set of all permutations, the algorithm has difficulties to converge for large $`N`$, while without any permutation it readily does. This dilemma can be solved by slowly interpolating from fully random permutations to the identity permutation: Execution of the main loop (steps 4-10) is dived into blocks of, e.g., 200 iterations. After terminating a block, the $`y_l(t)`$ are reset to the pair of series with the lowest $`\delta `$ so far, and the random permutations $`\mathrm{\Pi }_l`$ are replace by the permutations which rank order the sequences $`\{(t+h\nu _l(t))modN\}_t`$ ($`t=0,\mathrm{},N1`$). Here $`\{\nu _l(t)\}_t`$ are sequences of independent random numbers distributed equally between $`0`$ and $`1`$ which are renewed after each block and, as usual, $`umodN:=umN`$, where $`m`$ is the largest integer with $`umN0`$. The parameter $`h`$ controls the interpolation between fully random permutations ($`h\mathrm{}`$) and the identity ($`h=0`$). Successively smaller values of $`h`$ are used after each block. As an example for a pair of long time series ($`N=8192`$, $`k_1=k_2=24`$), Pair B was generated using $`h=400,200,100,0,0`$ (convergence is reached after 5 blocks) and all other parameters as for Pair A. For these long time series the fluctuations of the phases are clearly superimposed by a linear drift of the phase difference (Fig. 5). For *Pair A* this is not obvious, and one might be tempted to attribute the difference in frequency there to a “mere fluctuation”. A Matlab code that implements the algorithm in the form used to generate Pair B is available on the internet555http://www.fdm.uni-freiburg.de/$``$axel/onthelimits/.

## 4 The influence of end-point effects for two spectral estimators

Calculating the periodogram of a time series (and the empirical autocorrelation function as the Fourier transform thereof) is the correct first step toward estimating its power spectrum only if end-point effects can be neglected. This is the case when the time series allows a meaningful periodic continuation (which the algorithm presented here guarantees), or if the time series are sufficiently long. For two methods of spectral estimation it shall here be demonstrated that the condition of sufficient length can also be met by the generated pairs of time series – without making used of periodic continuation. That is, although the two methods differ in the way endpoints are treated, in the way the periodograms are averaged, and also slightly in their results, each of them gives virtually identical results for the two time series of a pair.

First consider Welch’s method (Welch, 1967). Both time series of Pair B are split into $`126`$ segments of $`192`$ samples length with $`2/3`$ overlap. The segments are tapered using a Hamming window and, in order to interpolate the resulting power spectra, extended to length 2048 by appending zeros. The periodograms of all segments of a time series are averaged. As is shown in Fig. 6, the estimated power spectra are nearly identical for both series. The differences are much smaller than then the difference of the oscillation frequencies.

Next, autoregressive spectral estimation is used. Applying Burg’s method (Marple, 1987), a $`p`$-th order autoregressive process

$`x(t)=\left({\displaystyle \underset{k=1}{\overset{p}{}}}a_kx(tk)\right)+e(t)`$ (8)
with uncorrelated, identically distributed $`e(t)`$, is fitted to the two time series of Pair B. The resulting estimates for the coefficients $`(a_1,\mathrm{},a_p)`$ for $`p=10`$ are listed in Table 1. No attempt was made to find the “optimal” model order $`p`$, since for all reasonable values of $`p`$ the same observation is made: The power spectra of the two fitted autoregressive processes are virtually indistinguishable. For example, with $`p=10`$ the spectral maxima are both located at the angular frequency $`\omega :=2\pi f=0.4887`$. They differ only by $`\mathrm{\Delta }\omega =2\times 10^6`$. In order to get an estimate for the errors of the autoregressive spectral estimates, the first and second half of the two time series are fitted separately. Then the enforced property of having equal periodograms is lost and the estimated spectra differ slightly. The maxima of the resulting power spectra are located at $`\omega =0.4876`$ and $`0.4897`$ for $`x_1(t)`$, and $`0.4890`$ and $`0.4884`$ for $`x_2(t)`$. For comparison, the oscillation frequencies of $`x_1(t)`$ and $`x_2(t)`$ are $`\omega =0.4725`$ and $`0.5093`$ respectively (see inset of Fig. 6). Hence, it is not possible to explain the undetectability of the difference in frequency between $`x_1(t)`$ and $`x_2(t)`$ of Pair B by a lack of spectral resolution.

## 5 Conclusion

Several examples were given that demonstrate that the power spectra of time series alone do not contain sufficient information to ensure that the oscillation frequencies (cycles per unit time) of the time series could precisely be determined. These examples illustrate a corresponding general theorem (Rossberg, 2002). Furthermore, it was shown that this result does not depend on the particular method of spectral estimation, as long as it is functional in the periodogram or the empirical autocorrelation function, respectively. There is also no other method of evaluating the periodogram which would yield the oscillation frequency. This is shown by introducing a general iterative method to generate pairs of time series of arbitrary length which have identical periodograms but different frequencies of oscillation.

## Figure Captions