text:
# Redshift estimation of clusters by wavelet decomposition of their Sunyaev-Zel’dovich morphology

## 1 Introduction

Inverse Compton scattering of cosmic microwave background (CMB) photons off thermal electrons within the hot intra-cluster medium (ICM) of galaxy clusters produce fluctuations in the surface brightness of the CMB, an effect known as the thermal Sunyaev-Zel’dovich (SZ) effect (e.g. Sunyaev & Zel’dovich, 1972, 1980; Rephaeli, 1995). Imaging clusters of galaxies through their SZ signature has, until recently, been a very challenging undertaking. To date, the development of detectors and new techniques have allowed high quality interferometric imaging of more than fifty clusters of galaxies (Carlstrom et al., 2002), despite incomplete coverage of the Fourier plane. In the foreseeable future, the availability of detectors in the microwave regime with angular resolutions surpassing $`1\mathrm{}`$ and sensitivities below $`\mu \mathrm{K}`$ (e.g., the South Pole Telescope, described in detail in Carlstrom et al., 2002), will probe the hot plasma in galaxy clusters out to large redshifts providing SZ based wide field galaxy cluster catalogues and yielding a multitude of information about cluster formationand the cosmological model (Birkinshaw, 1999).

In particular, the abundance of clusters as a function of redshift has been shown to be a very sensitive probe of the cosmological model (Eke et al., 1998; Henry, 2000). The near independence of the line-of-sight SZ amplitude on cluster redshift makes the SZ effect the main tool for detecting galaxy clusters at high redshifts, $`0.5\stackrel{<}{}z\stackrel{<}{}2`$ (the upper limit depends on cosmology quite sensitively). This range of redshifts is especially important for probing the nature of the dark energy of the universe, since during this era it is expected to evolve rapidly until it eventually dominates over the other cosmological fluids. In order to obtain precise constraints on cosmological models it is essential to have accurate measurements of the redshift distribution of galaxy clusters (see Haiman et al., 2001).

Normally, one determines the distance to the cluster by photometric or spectroscopic observations of the cluster member galaxies. Unfortunately, this is a very challenging and time consuming task, in particular, when one considers the very large number of mostly high redshift clusters expected to be observed with sensitive future SZ instruments – The Planck satellite alone is expected to detect about $`10^4`$ clusters (Bartelmann, 2001). In order to replace photometric follow-ups we aim at inferring the distance to a cluster from SZ data alone for a future generation of experiments with increased angular resolution of about $`1\mathrm{}`$.

Theoretically, the cold dark matter (CDM) hierarchical clustering paradigm predicts a universal profile for dark matter halos that depends only on two parameters: core radius and density (Navarro et al., 1995). In addition, the same theory provides a very simple recipe for the mass accretion history of a certain halo as a function of its formation and observation redshift (Wechsler et al., 2002; van den Bosch, 2002; Zhao et al., 2003). Using these relations together with simple assumptions like hydrostatic equilibrium and isothermality, one can expect that in the framework of the spherical collapse model the observable SZ flux and apparent size should provide measures of the cluster’s mass and distance.

Indeed, using scaling relations, Diego et al. (2003) have demonstrated the viability of determining reliable morphological redshifts and examined different SZ observables with respect to their distance sensitivity. Among those observables, they showed that the cluster apparent size and central amplitude are promising distance indicators, once their degeneracy is broken.

The main goal of this work is to derive redshifts of clusters based solely on their resolved SZ images by modeling the evolution of their structural parameters with redshift from the data set itself. This phenomenological approach does not depend on a priori assumptions about scaling relations that are valid only for spherically symmetric and relaxed systems.

Specifically, the structural morphology of the cluster’s pressure profile in an SZ observation is characterised by wavelet analysis.111There are also various ways of characterising the cluster’s density profile in an SZ observation that are more or less susceptible to noise, for instance the fitting of a $`\beta `$-profile (Cavaliere & Fusco-Femiano, 1978) to the electron density. We are able to show that there is a simple relation between the distribution of moments over various scales in wavelet space and the cluster properties which can be described with simple phenomenological functions. Furthermore, the parameters of these functions are shown to follow a well defined and simple redshift dependence. Wavelet analysis has been chosen because it maintains the scale and positional information of cluster morphology, hence, it makes isolation and suppression of various unwanted contributions to the observed signal possible while it reliabely upholds the underlying behavior. We note however, that Fourier space analysis could in principle yield very similar results.

Hydrodynamically simulated clusters are used to demonstrate the method and to set limits on the redshift uncertainty expected in this approach. The simulated clusters used in the analysis are close to virialisation, e.g. merging systems are excluded. Under this restriction, both the relation between the observed quantity and the cluster physical parameters as well as the structural parameters are well defined. In addition, simulated clusters ignore radiative and feedback processes, the effect of which is discussed later in the paper.

In the observational application, the evolution of the structural parameters following from wavelet decomposition could be calibrated from a (relatively small) learning set of high quality SZ clusters with known (photometric/spectroscopic) redshifts.

Our method relies crucially on the availability of resolved SZ cluster images. Therefore, throughout the paper we assume an instrumental resolution of $`1\mathrm{}`$, where massive clusters should be resolved even at the largest redshifts considered here. Indeed, future instruments such as the South Pole Telescope222http://astro.uchicago.edu/spt/ (Carlstrom et al., 2002) or the Atacama Cosmology Telescope333http://www.hep.upenn.edu/~angelica/act/act.html are designed to yield observations of up to $`10^4`$ galaxy clusters with masses $`\stackrel{>}{}10^{14}M_{}`$ ($`1\text{ }\mu \text{K}`$ sensitivity) and $`1\mathrm{}`$ resolution.

This article is organised as follows: After basic definitions concerning the SZ effect in Sect. 2 and wavelets in Sect. 3, the simulations are outlined in Sect. 4. The capability of wavelets with respect to distance estimation is examined in Sect. 5. Possible systematics are addressed in Sect. 6. A summary of the techniques in Sect. 7 and of the results in Sect. 8 concludes the article.

## 2 Sunyaev-Zel’dovich definitions

The SZ effect has been described in detail by many authors (for a comprehensive review see Birkinshaw, 1999); here we briefly review its main aspects. The SZ effect arises because CMB photons experience Compton-scattering off electrons of the dilute intra-cluster plasma. The CMB spectrum is modulated as photons are redistributed from the low-frequency part of the spectrum below $`218\text{ GHz}`$ to higher frequencies. The change in thermodynamic CMB temperature due to the thermal SZ effect is

$$\frac{\mathrm{\Delta }T}{T}(\mathit{\varphi })=y(\mathit{\varphi })\left(x\frac{e^x+1}{e^x1}4\right)2y(\mathit{\varphi })\text{for }x1\text{,}$$
(1)
where $`x=h\nu /kT_{\mathrm{CMB}}`$ is the dimensionless frequency. In the Rayleigh-Jeans limit ($`x1`$), the change in temperature is asymptotically equal to $`2y(\mathit{\varphi })`$. The SZ amplitude at location $`\mathit{\varphi }`$, which is known as the Comptonisation parameter $`y(\mathit{\varphi })`$, is defined as the line-of-sight integral of the temperature-weighted thermal electron density:

$$y(\mathit{\varphi })=\frac{\sigma _\mathrm{T}k}{m_\mathrm{e}c^2}dln_\mathrm{e}(\mathit{\varphi },l)T_\mathrm{e}(\mathit{\varphi },l)\text{.}$$
(2)
where $`m_\mathrm{e}`$, $`c`$ and $`k`$ denote electron mass, speed of light and Boltzmann’s constant, respectively. $`T_\mathrm{e}(\mathit{\varphi },l)`$ and $`n_\mathrm{e}(\mathit{\varphi },l)`$ are electron temperature and electron number density at position $`\mathit{\varphi }`$ and distance $`l`$.

## 3 Wavelets

### 3.1 Wavelet definitions

During the last decade, wavelet analysis has become a popular tool in various data analysis and image processing applications. The main appeal of wavelet functional bases stems from their simultaneous localisation of a signal in both the wavenumber and position domain, where they make orthogonal and complete projections on modes belonging to both spaces possible. In particular, the discrete wavelet families, by virtue of their constituting a complete basis, provide a unique and fast decomposition of the images into wavelet expansion coefficients. Statistics in terms of the $`q^{\mathrm{th}}`$ moments of the distribution of wavelet coefficients as a function of scale compresses the signal contained in an image into a small number of parameters and yields information surpassing that derived in traditional Fourier analysis.

Following Daubechies & Bates (1993) and Muzy et al. (1993), the wavelet transform of a 2-dimensional image is defined as a convolution of the function $`y(𝒙)`$ to be analysed with the wavelet $`\psi _\sigma \left(\left|𝒙𝝁\right|\right)`$:

$$\chi (𝝁,\sigma )=\mathrm{d}^2xy(𝒙)\psi _\sigma \left(\left|𝒙𝝁\right|\right)\text{.}$$
(3)
High values for $`\chi (𝝁,\sigma )`$ are obtained in case of a match between the features of $`y(𝒙)`$ and the wavelet $`\psi _\sigma (𝒙)`$ at position $`𝝁`$ and scale $`\sigma `$. From the wavelet expansion coefficients $`\chi (𝝁,\sigma )`$ on scale $`\sigma `$ at location $`𝝁`$ one obtains the wavelet moments $`X_q(\sigma )`$ by integration over all displacements $`𝝁`$:

$$X_q(\sigma )=\mathrm{d}^2\mu \left|\chi (𝝁,\sigma )\right|^q\text{.}$$
(4)
The exponent $`q`$ defines the order of the wavelet moment $`X_q(\sigma )`$. Values for $`q`$ equal or larger than 2 allow noise suppression. The logarithm $`\mathrm{ln}X(\sigma ,q)`$ of the wavelet moment as a function of logarithmic scale $`\mathrm{ln}\sigma `$ constitutes the wavelet spectrum. The $`X_q(\sigma )`$-statistic is the main tool used in this study for characterising the morphology of SZ clusters.

### 3.2 Application of wavelets to a cluster profile

#### 3.2.1 Analytic wavelet transform of a cluster $`y`$-profile

In order to illustrate our idea of determining cluster sizes via wavelet decomposition, the wavelet transform of a King profile, which is known to describe the SZ morphology of clusters to first order, is performed. As an analysing wavelet, the Mexican-hat wavelet was chosen for simplicity.

It is favourable to compute the convolution in the definition of $`\chi (𝝁,\sigma )`$ in the Fourier domain. By virtue of eqn. (6),

$`\chi (𝝁,\sigma )`$ $`=`$ $`{\displaystyle \mathrm{d}^2xy(𝒙)\psi _\sigma \left(𝒙𝝁\right)}`$ (5)
$`=`$ $`(2\pi )^2{\displaystyle \mathrm{d}^2kY(𝒌)\mathrm{\Psi }_\sigma (𝒌)\mathrm{exp}(\mathrm{i}𝒌𝝁)\text{,}}`$ (6)
the convolution reduces to a mere multiplication of the Fourier transforms $`Y(𝒌)`$ and $`\mathrm{\Psi }_\sigma (𝒌)`$ of the image $`y(𝒙)`$ and the wavelet $`\psi _\sigma (𝒙)`$, respectively. Restricting the order of the wavelet moment to $`q=2`$ and inserting the convolution theorem (6) into the definition (4) yields:

$`X_2(\sigma )`$ $`=`$ $`(2\pi )^4{\displaystyle \mathrm{d}^2\mu \left|\mathrm{d}^2kY(𝒌)\mathrm{\Psi }_\sigma (𝒌)\mathrm{exp}(\mathrm{i}𝒌𝝁)\right|^2}`$ (7)
$`=`$ $`(2\pi )^6{\displaystyle \mathrm{d}^2k\left|Y(𝒌)\right|^2\left|\mathrm{\Psi }_\sigma (𝒌)\right|^2\text{,}}`$ (8)
where the replacement $`\left|\mathrm{\Psi }_\sigma (𝒌)\right|^2=\left|\mathrm{\Psi }_\sigma (𝒌)\right|^2`$ holds for real wavelets.

The Mexican-hat wavelet is defined as the negative Laplacian of a Gaussian:

$$\psi _{\mathrm{MH}}(𝒙)=\psi _{\mathrm{MH}}(x)=_𝒙^2\left[\frac{1}{2\pi \sigma ^2}\mathrm{exp}\left(\frac{𝒙^2}{2\sigma ^2}\right)\right]\text{,}$$
(9)
whereof the Fourier transform $`\mathrm{\Psi }_{\mathrm{MH}}(𝒌)`$ is derived by twofold partial integration:

$`\mathrm{\Psi }_{\mathrm{MH}}(k)`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}^2x}{(2\pi )^2}\psi _{\mathrm{MH}}(𝒙)\mathrm{exp}(\mathrm{i}𝒌𝒙)}`$ (10)
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^2\sigma ^6}}{\displaystyle rdr(2\sigma ^2r^2)\mathrm{exp}\left(\frac{r^2}{2\sigma ^2}\right)J_0(kr)}`$ (11)
$`=`$ $`{\displaystyle \frac{k^2}{(2\pi )^2}}\mathrm{exp}\left({\displaystyle \frac{k^2\sigma ^2}{2}}\right)\text{,}`$ (12)
where the azimuthal symmetry and the definition of the zeroth order Bessel function of the first kind, $`2\pi J_0(kr)=_0^{2\pi }d\varphi \mathrm{exp}\left(\mathrm{i}kr\mathrm{cos}\varphi \right)`$ was used in the first step. Thus, the Fourier transform of the wavelet, $`\mathrm{\Psi }_\sigma (k)`$, is given by the Hankel transform of the Laplacian of a Gaussian.

For the determination of $`Y(𝒌)`$, we assume that the projected thermal electron density can be described by a spherically symmetric King profile, i.e. a $`\beta `$-model (Cavaliere & Fusco-Femiano, 1978) with $`\beta =1`$, core radius $`r_c`$ and central value of the Comptonisation parameter $`y_0`$:

$$y(𝒙)=y(r)=y_0\left[1+\left(\frac{r}{r_c}\right)^2\right]^1\text{.}$$
(13)
Then, the Fourier transform is easily computed:

$`Y(k)`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}^2x}{(2\pi )^2}y(𝒙)\mathrm{exp}(\mathrm{i}𝒌𝒙)}`$ (14)
$`=`$ $`{\displaystyle \frac{y_0r_c^2}{2\pi }}{\displaystyle dr\frac{r}{r_c^2+r^2}J_0(kr)}={\displaystyle \frac{y_0r_c^2}{2\pi }}K_0(kr_c)\text{,}`$ (15)
where in eqn. (15) the definition of the zeroth order modified Bessel function of the second kind $`K_0(kr_c)`$ was inserted.

Substituting eqns. (12) and (15) into eqn. (8) and exploiting the azimuthal symmetry of the functions $`y(𝒙)`$ and $`\psi (𝒙)`$ yields an analytic integral for $`X_2(\sigma )`$:

$$X_2(\sigma )=2\pi y_0^2r_c^4_0^{\mathrm{}}dkk^5\mathrm{exp}\left(\sigma ^2k^2\right)K_0^2(kr_c)\text{.}$$
(16)
After evaluation of the integral in eqn. (16), the wavelet transform of the $`\beta `$-profile reads as follows:

$`X_2(\sigma )`$ $`=`$ $`{\displaystyle \frac{\pi ^{3/2}y_0^2}{2r_c^2}}\alpha ^6𝒢_{2,3}^{3,1}\left(\alpha ^2|\genfrac{}{}{0pt}{}{2\frac{1}{2}}{\mathrm{0\hspace{0.17em}0\hspace{0.17em}0}}\right)\text{,}`$ (17)
where $`\alpha =r_c/\sigma `$ has been substituted. The function $`𝒢`$ is Meijer’s G-function, the exact definition of which is given by Gradshteyn & Ryzhik (1994). It is an interesting consistency to note that apart from the normalisation, the functional shape of eqn. (17) only depends on $`\alpha `$, i.e. on the core radius $`r_c`$ expressed in units of the wavelet scale $`\sigma `$.

#### 3.2.2 Asymptotics of the analytical wavelet transform

The asymptotic behavior of $`X_2(\sigma )`$ at the limit of $`\sigma r_c`$ can be explored by substituting the expressions given in eqns. (13) and (9) into eqn. (6), and exchanging, by partial integration, the function on which the Laplacian operates. In the limit of interest the Gaussian can be replaced by a Dirac-$`\delta `$ distribution. Substituting all of this into eqn. (4) yields that $`lim_{\sigma 0}X_2(\sigma )`$ is proportional to $`y_0^2`$ and independent of $`\sigma `$, i.e. the normalisation of the wavelet spectrum measures the square of the central Comptonisation parameter $`y_0`$:

$$X_2(\sigma )=\frac{32\pi }{15}\frac{y_0^2}{r_c^2}\text{for}\sigma r_c\text{.}$$
(18)
In the opposite limit, i.e. $`r_c\sigma `$, one can use the fact that the King-profile is highly peaked at the center and that it is convolved with a Mexican-hat wavelet guaranteeing the convergence of the integral in eqn. (16) at $`\mathrm{}`$. In the limit of $`r_c0`$ this integral is dominated by the value at $`k=0`$. Therefore, one can approximate the King-profile with a Dirac-$`\delta `$ distribution and show that asymptotically the $`lim_{r_c0}X_2(\sigma )`$ is proportional to $`\sigma ^6`$:

$$X_2(\sigma )\frac{y_0^2r_c^4}{\sigma ^6}\text{for}\sigma r_c\text{.}$$
(19)
The sensitivity of the wavelet spectrum $`X_2(\sigma )`$ on cluster size $`r_c`$ is illustrated by Fig. 6. The wavelet spectrum is constant for $`\sigma r_c`$, has an $`r_c`$-dependent break and drops off asymptotically $`\sigma ^6`$ for $`\sigma r_c`$. Naturally, the scale $`\sigma `$ at which the transition from one asymptotic regime to the other occurs, is determined by the value of $`r_c`$, i.e. the cluster size.

Motivated by this example, the wavelet moments $`X_q(\sigma )`$ obtained from real data (Sect. 5.3) will be fitted with a power law with an exponential cutoff, where the cutoff indicates the cluster size and the amplitude is proportional to some power of the central Comptonisation parameter $`y_0`$.

#### 3.2.3 Finite instrumental resolution

The influence of finite instrumental resolution can easily be incorporated by an additional factor $`\left|B(k)\right|^2`$ in eqn. (8):

$$X_2(\sigma )=(2\pi )^6\mathrm{d}^2k\left|Y(k)\right|^2\left|\mathrm{\Psi }_\sigma (k)\right|^2\left|B(k)\right|^2\text{,}$$
(20)
where $`B(k)`$ is the Fourier transform of the (azimuthally symmetric) beam profile $`b(x)`$, which is for simplicity assumed to be of Gaussian shape with FWHM $`=\sqrt{8\mathrm{ln}(2)}\sigma _b`$:

$`B(k)`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}^2x}{(2\pi )^2}b(x)\mathrm{exp}(\mathrm{i}𝒌𝒙)\text{ with}}`$ (21)
$`b(𝒙)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \sigma _b^2}}\mathrm{exp}\left({\displaystyle \frac{𝒙^2}{2\sigma _b^2}}\right)\text{.}`$ (22)
This effectively replaces $`\sigma `$ in eqn. (17) by the harmonic mean $`\sqrt{\sigma ^2+\sigma _b^2}`$, which limits the range of accessible wavelet scales to $`\sigma >\sigma _b`$.

### 3.3 Analogy to power spectra in Fourier analysis

By interpreting the wavelet spectrum in eqn. (8) as the variance of the fluctuations on the scale $`\sigma `$, one may draw an analogy to Fourier decomposition:

$$\mathrm{var}\left[y(x)\right]=X_2(\sigma )=(2\pi )^4\mathrm{d}^2kP(k)\left|\mathrm{\Psi }_\sigma (k)\right|^2\text{,}$$
(23)
where $`P(k)=(2\pi )^2\left|Y(k)\right|^2`$ is the Fourier power spectrum. The wavelet $`\psi (x)`$ now adopts the role of a filter function on scale $`\sigma `$. This filter function reads in real space in the case of the Mexican hat wavelet:

$$\psi _{\mathrm{MH}}(𝒙)=\frac{2\sigma ^2𝒙^2}{2\pi \sigma ^6}\mathrm{exp}\left(\frac{𝒙^2}{2\sigma ^2}\right)\text{.}$$
(24)
Therefore, our method is equivalent to considering power spectral analysis of filtered fields and higher order Fourier space moments.

## 4 Simulations

The accuracy in the determination of redshift $`z`$ was assessed by examining the performance on numerical simulations: First, simulations of cosmological structure formation including gas physics have been carried out in order to model the evolution of clusters (Sect. 4.1). Subsequently, maps of the Compton-$`y`$ parameter have been produced by using an interpolation kernel with an adaptive smoothing length for projecting the Compton-$`y`$ parameter along the line-of-sight (Sect. 4.2). By applying selection criteria favouring virialised systems a cluster sample was compiled (Sect. 4.3). Finally, aiming at realistic single frequency SZ observations, we simulated the ambient CMB fluctuations that act as the primary source of noise (Sect. 4.4) and combined the resulting realisations of the CMB with the cluster maps (Sect. 4.5).

The assumed cosmological model is the standard $`\mathrm{\Lambda }`$CDM cosmology, which has recently been supported by the WMAP satellite (Bennett et al., 2003; Spergel et al., 2003). Parameter values have been chosen as $`\mathrm{\Omega }_\mathrm{M}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`H_0=100h\text{km}\text{s}^1\text{ Mpc}^1`$ with $`h=0.7`$, $`\mathrm{\Omega }_\mathrm{B}=0.04`$, $`n_\mathrm{s}=1`$ and $`\sigma _8=0.9`$.

### 4.1 SPH cluster simulations

A simulation of cosmological structure formation kindly provided by V. Springel and L. Hernquist (Springel & Hernquist, 2002; White et al., 2002) constitutes the basis of our analysis. In a cubical box of comoving side length $`100\text{ Mpc}/h`$ with periodic boundary conditions a smoothed particle hydrodynamic (SPH) simulation comprising $`216^3`$ dark matter particles as well as $`216^3`$ gas particles was run and snapshots were saved at 23 redshifts ranging from $`z=0.102`$ out to $`z=1.114`$. The comoving spacing along the line-of-sight of two subsequent boxes is $`100\text{ Mpc}/h`$. Purely adiabatic gas physics and shock heating were included, but radiative cooling and star formation were ignored, which however does not result in significant differences in SZ morphology, as shown by White et al. (2002) but impacts on the scaling relations as demonstrated by da Silva et al. (2001).

Overdensities are identified using a friends-of-friends algorithm with the linking length $`b=0.164`$, which yields all member particles of a cluster in conjunction with a spherical overdensity code, from which virial quantities are estimated. We computed the mass $`M_{\mathrm{vir}}`$ inside a sphere of radius $`r_{\mathrm{vir}}`$, interior to which the average density was 200 times the critical density $`\rho _{\mathrm{crit}}=3H_0^2/(8\pi G)`$. The angle subtended by twice the virial radius is denoted as $`\theta _{\mathrm{vir}}`$. We imposed a lower mass threshold of $`M_{\mathrm{vir}}510^{13}M_{\mathrm{}}/h`$.

### 4.2 SZ map preparation

Square maps of the Compton-$`y`$ parameter of the selected clusters were generated by SPH projection of all member gas particles onto a cubical grid with $`128^2`$ mesh points. The (comoving) side length $`s`$ of the map was adapted to the cluster size, such that the comoving resolution $`g=s/128`$ of the grid is specific to a given map. Examples of Sunyaev-Zel’dovich maps are given in Fig. 2.

If the particle $`p`$ at position $`𝒓_p=(x_p,y_p,z_p)`$ has a smoothing length $`h_p`$, an SPH electron number density estimate $`n_p`$, and an SPH electron temperature $`T_p`$, the Compton-$`y`$ parameter at the pixel at position $`𝒙`$ is given by:

$$y(𝒙)=\frac{\sigma _\mathrm{T}k}{m_ec^2}\frac{h_p^3}{g^2}\underset{p}{}\left[\underset{xg/2}{\overset{x+g/2}{}}dx_p\underset{yg/2}{\overset{y+g/2}{}}dy_p\underset{h_p}{\overset{h_p}{}}dz_pw\left(\frac{r}{h_p}\right)n_pT_p\right]$$
(25)
$$\text{with }r=\sqrt{(x_px)^2+(y_py)^2+z_p^2}\text{.}$$
(26)
Here, we assumed complete ionisation and primordial element composition of the ICM for the determination of electron number density and temperature. In this way we produced projections along each of the three coordinate axes. The function $`w`$ is the spherically symmetric cubic spline kernel suggested by Monaghan & Lattanzio (1985), which was also used in the SPH simulation:

$$w(u)=\frac{8}{\pi }\{\begin{array}{cc}16u^2+6u^3,\hfill & 0u1/2\hfill \\ 2(1u)^3,\hfill & 1/2<u1\hfill \\ 0,\hfill & u>1\hfill \end{array}\text{ with }u=r/h_p\text{.}$$
(27)
The fact that the kernel $`w`$ is defined on a compact support $`u\left[0\mathrm{}1\right]`$ greatly reduces the computational effort.

### 4.3 Cluster selection

Clearly, the wavelet redshift estimation relies on the clusters not being in the state of violent merging. Apart from the minimal mass of $`M_{\mathrm{min}}=510^{13}M_{\mathrm{}}/h`$, that translates into a minimally required line-of-sight Comptonisation amplitude $`y_{\mathrm{min}}`$, clusters have been selected in order to meet the following prerequisites:

* By visual inspection it was made sure that the clusters show but a single peak in the Compton-$`y`$ map in order to exclude systems in the late phase of a merger.
* The SZ morphology is required not to be too elongated. By fitting a 2-dimensional $`\beta `$-model $`y_\beta (𝒙)`$ to the SZ profile $`y_{\mathrm{data}}(𝒙)`$, values for the semi-axes $`r_x`$ and $`r_y`$ are derived. $`90\%`$ of the clusters within the selected sample have axis ratios $`q=r_y/r_x`$ smaller than 0.8 and ellipticities $`e=\sqrt{r_x^2r_y^2}/r_x`$ below 0.6.
* Residual deviations from the canonical $`\beta `$-profile ought to be small. The $`rms`$-deviation $`v`$ of the cluster from the best-fitting $`\beta `$-profile,

$$v=\sqrt{\left(\frac{y_{\mathrm{data}}(𝒙)y_\beta (𝒙)}{y_\beta (𝒙)}\right)^2_𝒙}\text{,}$$
(28)
was smaller than $`25\%`$ for $`90\%`$ of our clusters.

Applying these selection criteria, 10 clusters were selected from each of the 23 redshift bins, yielding with the three orthogonal projections of each cluster a total number of 690 maps with which the accuracy of the wavelet method in estimating redshifts was assessed. The distributions of the ellipticities $`e`$ and the integrated residuals $`v`$ are shown in Fig. 3. The same distributions were derived for the smoothed cluster maps, were the effects of finite instrumental resolution have been incorporated. As Fig. 3 suggests, the beam does not have a major impact on the morphological properties of the cluster sample, which is due to its narrowness of only $`1\mathrm{}`$ (FWHM).

### 4.4 CMB map generation

CMB anisotropies are assumed to be a particular realisation of a Gaussian random field. Aiming at simulating a realisation of the CMB on a square, flat map, we take temperature fluctuations $`\theta (\mathit{\varphi })`$ relative to the average CMB temperature of $`T=2.726\mathrm{K}`$ to be the independent random field,

$$\theta (\mathit{\varphi })\frac{T(\mathit{\varphi })T}{T}\text{.}$$
(29)
The flat, two-dimensional power spectrum $`P_\theta (\mathrm{})`$ is defined via:

$$\mathrm{\Theta }(\mathbf{})\mathrm{\Theta }^{}(\mathbf{}^{})(2\pi )^2\delta _D(\mathbf{}\mathbf{}^{})P_\theta (|\mathbf{}|)\text{,}$$
(30)
where $`\mathrm{\Theta }(\mathbf{})`$ denotes the Fourier transform of $`\theta (\mathit{\varphi })`$. The simulation of the CMB temperature fluctuations on a flat square map now consists of the following two steps:

* The angular power spectrum $`C_{\mathrm{}}`$ is computed for the flat $`\mathrm{\Lambda }`$CDM-Universe using the CMBfast code by Seljak & Zaldarriaga (1996). In addition to the cosmological parameters being already described in Sect. 4, we use adiabatic initial conditions and set the primordial He-mass fraction to $`X_{\mathrm{He}}=0.24`$ and the Thomson optical depth to $`\tau =0.17`$ (Spergel et al., 2003). The angular power spectrum of the CMB is normalised to COBE data. Since the SZ effect distorts the CMB only on small angular scales, the flat sky approximation $`\mathrm{}1`$ is fulfilled and it is appropriate to replace the spherical harmonics with plane waves. Hu (2000) has shown that the 2-dimensional flat power spectrum $`P_\theta (\mathrm{})`$ is approximately equal to its angular analogue: $`C_{\mathrm{}}P_\theta (\mathrm{})`$.
* Then, Gaussian random variables are generated on a complex two-dimensional grid in Fourier space with variance $`\sigma ^2(\mathrm{})=P_\theta (\mathrm{})`$ according to the absolute value of their wave vectors $`\mathbf{}`$. Inverse Fourier transform yields a realisation of the temperature anisotropies $`\theta (\mathit{\varphi })`$.

### 4.5 Simulated single-frequency SZ observations

For SZ clusters observed with a single-frequency interferometer (e.g., the CBI experiment, Halverson et al., 2002)111http://www.astro.caltech.edu/~tjp/CBI/, it is important to examine the applicability of the $`X_q(\sigma )`$-statistic. For the purpose of this paper, it suffices to consider observations at small frequencies $`\nu `$. Thus, the Compton-$`y`$ maps are combined with realisations of the CMB fluctuations by taking advantage of eqn. (1) in the Rayleigh-Jeans limit,

$$T(\mathit{\varphi })=\left[12y(\mathit{\varphi })\right]\left[1+\theta (\mathit{\varphi })\right]T\text{.}$$
(31)
Fig. 4 shows the Compton-$`y`$ map of a nearby cluster of $`2.210^{14}M_{\mathrm{}}/h`$ at redshift $`z=0.102`$ combined with a patch of the CMB constructed by the algorithm described above. In this map, the average CMB temperature $`T`$ was subtracted. In order tomimic observations, the resulting combined maps are smoothed with a Gaussian beam with FWHM of $`\sqrt{8\mathrm{ln}(2)}\sigma _b=1\mathrm{}`$.

In the case of multi-frequency SZ observations the SZ signature can be easily distinguished from the CMB signal. Therefore, for these cases the CMB background is ignored and not included in the simulated cluster SZ images. Nevertheless, finite instrumental resolution was taken care of and the SZ maps were convolved with a Gaussian kernel of $`\sqrt{8\mathrm{ln}(2)}\sigma _b=1\mathrm{}`$ (FWHM).

## 5 Analysis

In this section, the analysis is explained step by step: After introducing the wavelet families (Sect. 5.1), the wavelet spectrum and the parameters deduced from it are described (Sect. 5.2 and Sect. 5.3). The correlations of the wavelet spectral parameters with physical quantities are discussed (Sect. 5.4). The measurement principle and the breaking of degeneracy is illustrated in Sect. 5.5. Successively, the intercorrelation of the wavelet parameters and the shape of the parameter space is explored by principal component analysis (Sect. 5.6). Then, gauge functions for modelling the redshift dependence ofthe parameters are proposed (Sect. 5.7). An important issue for single-frequency interferometers is the influence of primary CMB fluctuations on the wavelet spectrum and their suppression (Sect. 5.8). Finally, the redshift of the clusters are estimated by maximum likelihood techniques (Sect. 5.9).

### 5.1 Wavelet basis functions

In the analysis a wide range of wavelets with different functional shapes was employed, although the symlet wavelet basis introduced by Daubechies & Bates (1993) yielded particularly good results. Due to their symmetry and peakiness, symlets are seemingly especially suited for analysing SZ morphologies, because spiky features are useful for edge detection, and since apparent cluster size is the primary observable, the usage of a peaked wavelet is justified. Other wavelet families that found application in our analysis were Daubechies’ wavelets, coiflets and biorthogonal wavelets. Fig. 5 compares the functional shape of the different wavelet families.

The analysis proceeds by measuring the wavelet moments on smoothed comoving maps of the Compton-$`y`$ parameter following the definition in Sect. 3. The scale $`\sigma `$ of the resulting wavelet spectrum is then converted to angular units. Because our SZ maps are computed on a grid of $`128^2`$ mesh points with adaptively chosen side length for each cluster, our dynamical range of the wavelet spectra always comprises approximately two decades. However, this is no fundamental limitation of this approach because the maps are featureless below the smoothing scale of $`1\mathrm{}`$ (FWHM).

### 5.2 Measurement of wavelet quantities

In order to derive the actual flux decrement or, equivalently, the decrement in antenna temperature from the line-of-sight Compton-$`y`$ amplitude, its value at each pixel needs to be multiplied with the solid angle it subtends. For the conversion, a standard $`\mathrm{\Lambda }`$CDM-cosmology was assumed, the parameters of which have already been given in Sect. 4. Thus, the pixel amplitudes were modified according to:

$$y(𝒙)y(\mathit{\varphi })=y(𝒙)4\mathrm{arctan}^2\left[\frac{g}{2w(z)}\right]\text{,}$$
(32)
where $`w(z)`$ is the comoving distance in the model cosmology and $`g`$ denotes the comoving size of a single pixel. It should be emphasised that the wavelet coefficients $`\chi (𝝁,\sigma )`$ are evaluated on a comoving grid, which has been adapted to the cluster size before converting the wavelet scale $`\sigma `$ to angular units. This, however, should not pose a problem for real observations, provided the sampling scale is of the same order of magnitude as the angular scale of our finest pixels.

In order to obtain dimensionless quantities, the unit of the wavelet $`\psi _\sigma (𝒙)`$ has been set to inverse steradians, such that the wavelet expansion coefficients $`\chi (𝝁,\sigma )`$ and the wavelet moments $`X_q(\sigma )`$ are dimensionless, irrespective of $`q`$. For numerical convenience, the pixel amplitudes in the combined SZ maps have been multiplied with $`10^{12}`$.

The summation in the definition of the wavelet moment $`X_q(\sigma )`$ in eqn. (4) discards the information about the position $`𝝁`$ at which the wavelet expansion coefficient $`\chi (𝝁,\sigma )`$ is evaluated. Consequently, the position of a cluster inside the observing frame does not influence the wavelet decomposition.

### 5.3 Wavelet spectrum of SZ cluster maps

Due to the lack of any analytical generalisation of eqn. (17) for $`q2`$, deviations of the Compton-$`y`$ map from a King profile and wavelets other than the simple Mexican hat, we decided to explore phenomenological functions for describing the wavelet spectrum. The simplicity of the shape of the wavelet spectrum shown in Fig. 6 implies that the model function,

$$\mathrm{ln}X_q(\sigma )a+s\mathrm{ln}\left(\sigma /\sigma _0\right)\sigma /c\text{,}$$
(33)
is able to extract all apparently contained information, i.e. the spectrum is described by means of three quantities: the amplitude $`a`$, the slope $`s`$ and a break at $`c`$. The parameter $`\sigma _0`$ has been included in eqn. (33) in order to obtain a formula that is dimensionally correct, although it does not yield any new information and this specific degree of freedom is already described by the variable $`a`$.

The usage of eqn. (33) implicitly discards information about asphericity and effectively determines an average of the cluster’s extension along its major axes. The problem would be significantly complicated by including asymmetry and considering vectorial nature of $`\sigma `$ (see Zaroubi et al., 1998, 2001).

Because the cutoff parameter $`c`$ is of great importance to our analysis, it needs to be derived reliably. Thus, the order of wavelet moments $`q`$ was restricted to $`q3`$, because larger $`q`$-values facilitate the determination of $`c`$. From Fig. 6 it is obvious that an increase in $`q`$ suppresses the value of $`X_q(\sigma )`$ at small scales $`\sigma `$ such that the curve develops a maximum in the vicinity of $`c`$. Additionally, by the choice of large values for $`q`$, the wavelet expansion coefficients $`\chi (\mu ,\sigma )`$ dominated by CMB noise are suppressed relative to those obtained in the central part of the cluster and consequently higher order wavelet moments $`X_q(\sigma )`$ provide a cleaner measurement. The range of sensible $`q`$-values is restricted by the fact that for increasing $`q`$ the moment $`X_q(\sigma )`$ is successively dominated by the largest wavelet expansion coefficient $`\chi (\mu ,\sigma )`$ and does no longer contain information of the structure to be analysed. In order to stabilise the fitting procedure we interpolate in between the wavelet moments $`X_q(\sigma )`$. This is justified because we expect a smooth variation of the wavelet spectrum according to Sect. 3.2.1.

### 5.4 Correlations with physical quantities

The parameters derived from the fit to the spectrum of wavelet coefficients have a physical interpretation. As illustrated in Sect. 3.2.1, the wavelet spectrum breaks at the cluster scale.

Therefore, one expects a correlation between the angular size of the cluster $`\theta _{\mathrm{vir}}`$ and the cutoff $`c`$, as shown by Fig. 7. Increasing weighting exponents $`q`$ shift the regression line to smaller values of $`c`$, which can be understood by the fact that larger values of $`q`$ suppress small wavelet expansion coefficients arising at the outskirts of the cluster, which in turn leads to a break in the wavelet spectrum at smaller scales.

Similarly, the amplitude $`a`$ determined by the fit is proportional to the integrated Compton-$`y`$ flux,

$$\mathrm{{\rm Y}}=\mathrm{d}^2\varphi y(\mathit{\varphi })=\frac{kT_{\mathrm{vir}}}{m_ec^2}\frac{\sigma _T}{d_\mathrm{A}(z)^2}\frac{1+f_H}{2}f_B\frac{M_{\mathrm{vir}}}{m_\mathrm{p}}\text{,}$$
(34)
as illustrated by Fig. 8. Here, $`f_b`$ denotes the baryon fraction, $`f_H`$ the hydrogen fraction, which determines the elemental composition and has been set to the primordial value of 0.76, and $`m_\mathrm{p}`$ is the proton mass. $`d_\mathrm{A}(z)`$ is the angular diameter distance in our cosmology.

The normalisation $`a`$ of the wavelet moments $`X_q(\sigma )`$ shows a steeper dependence on the integrated Comptonisation parameter $`\mathrm{{\rm Y}}`$ for larger choices of $`q`$, which is explained by the following argument: The amplitude $`a(q)`$ reflects the normalisation of the wavelet moments $`X_q(\sigma )`$. The integral in eqn. (4) is dominated by the largest wavelet expansion coefficient $`\chi (𝝁,\sigma )`$, taken to the $`q^{\mathrm{th}}`$ power. On the other hand, the wavelet expansion coefficients $`\chi (𝝁,\sigma )`$ are proportional to integrated Comptonisation $`\mathrm{{\rm Y}}`$, resulting in observed relation $`\mathrm{ln}\left[X_q(\sigma )\right]aq\mathrm{ln}(\mathrm{{\rm Y}})`$.

The slope $`s`$ is a measure of instrumental smoothing: Placing the same cluster at different redshifts would result in a blurred image of the more distant one. Keeping in mind that there is a close analogy between wavelet- and Fourier-transforms (as explained in Sect. 3.3), the wavelet moment $`X_q(\sigma )`$ as a function of $`\sigma `$ can be interpreted as the variance of the wavelet-filtered field. The instrumental beam introduces an additional filtering to the Compton-$`y`$ map (compare Sect. 3.2.3) and would cause the Fourier spectrum to drop at smaller values of the wave vector $`k`$, because the instrumental beam constitutes effectively a low-pass filter that is erasing structures smaller than its extension. Comparing clusters at different redshifts, it is clear that the drop in power happens at smaller scales in the case of the more distant cluster. Then the slope $`s`$, defined as $`\mathrm{d}\mathrm{ln}X_q(\sigma )/\mathrm{d}\mathrm{ln}\sigma `$ for $`\sigma r_c`$, is larger in the case of a resolved cluster compared to an unresolved cluster. This measure of the influence of finite instrumental smoothing varies only by a factor of two in the redshift and mass range considered here, but nevertheless serves as an indicator of cluster distance.

### 5.5 Measurement principle

Now, it is necessary to illustrate how a measurement of the total Comptonisation $`\mathrm{{\rm Y}}`$ and of the angular size $`\theta _{\mathrm{vir}}`$ suffices to derive a distance estimate. For that purpose, clusters are placed at unit distance and the distance dependences of the wavelet amplitude $`a`$ and the the cutoff $`c`$ are removed by the following formulae, since $`a`$ is a logarithmic measure of flux inside an solid angle element $`\mathrm{{\rm Y}}`$ and $`c`$ is a logarithmic measure of of angular extension $`\theta _{\mathrm{vir}}`$:

$`a_0`$ $`=`$ $`a(z)+2\mathrm{ln}\left(d_a(z)\right)`$ (35)
$`c_0`$ $`=`$ $`c(z)+\mathrm{ln}\left(d_a(z)\right)\text{.}`$ (36)
By applying simple scaling arguments, one expects the ratio $`a_0/c_0`$ to be equal to 5: From the wavelet amplitude $`a`$ one obtains $`a_0\mathrm{ln}\left(\mathrm{{\rm Y}}d_A(z)^2\right)\mathrm{ln}\left(M_{\mathrm{vir}}T_{\mathrm{vir}}\right)`$. Furthermore, from the spherical collapse model follows, that $`T_{\mathrm{vir}}M_{\mathrm{vir}}^{2/3}`$ (Navarro et al., 1995), which yields, together with $`M_{\mathrm{vir}}r_{\mathrm{vir}}^3`$, the relation $`a_0\mathrm{ln}\left(r_{\mathrm{vir}}^5\right)`$. Substituting $`c_0\mathrm{ln}\left(r_{\mathrm{vir}}\right)`$ gives the final result $`a_0/c_0=5`$.

Fig. 9 nicely illustrates how the degeneracy is broken and how a simple measurement of flux and angular extension suffices in order to derive a distance estimate: A crude fit to the distance corrected wavelet amplitude $`a_0`$ as a function of distance corrected wavelet cutoff parameter $`c_0`$ yields a slope of approximately 5.8, which corresponds well to the slope of $`5`$ expected from the theoretical consideration outlined above. If, hypothetically, the ratio $`a_0/c_0`$ was equal to 2, the measurements of flux and angular size would be completely degenerate and would not yield any distance information. This case corresponds to disks of equal surface brightness, where measurements of flux and angular size are completely degenerate and do not yield any distance information at all. It should be noted, that by adopting the usual scaling relations, one introduces a systematic error in slope that can amount to $`20\text{\%}`$ error.

### 5.6 Principal component analysis

In order to rate the extent to which the parameters derived in the fit to the spectrum of wavelet moments are independent, a principal component analysis (PCA, see e.g. Deeming, 1964) was performed. The PCA is determining a transformation to a new orthogonal coordinate system in parameter space spanned by $`a`$, $`c`$ and $`s`$, such that the variance along the first axis is maximised.

The first eigenvector of the matrix that describes the change of basis by the PCA reads $`x_{\mathrm{PCA}}=(0.65,0.70,0.32)`$ which has been derived for the spectral parameters for $`q=3`$ and with the sym2-wavelet as the analysing wavelet. The values similar in magnitude suggest that the variation in the data set is contained in all three parameters $`a`$, $`c`$ and $`s`$.

As can be read off from table 1, the parameter space is tightly constrained and all three parameters are interrelated, such that the data points form a narrow ray in parameter space. This result holds irrespective of the choice of $`q`$, although the scatter increases with higher choices for $`q`$. The values in table 1 have been determined without considering CMB fluctuations. Given the physical interpretation of the wavelet amplitude $`a`$ and the cutoff $`c`$, it is obvious that the tight correlation can be traced back to the self-similarity of clusters and the cluster scaling relations linking $`T_{\mathrm{vir}}`$, $`M_{\mathrm{vir}}`$ and $`r_{\mathrm{vir}}`$ that follow from the spherical collapse model. The scaling relations for SZ quantities derived by da Silva et al. (2003) support this view. This shows together with Sect. 5.4 and Sect. 5.5 hat both the cutoff $`c`$ and the amplitude $`a`$ are functions of $`M_{\mathrm{vir}}`$. By only considering systems close to virialisation we thus expect a fundamental narrow ray in parameter space which can be described by a single principal component.

### 5.7 Redshift dependence of the wavelet parameters

The parameters $`a`$ and $`c`$ are expected to decrease with increasing redshift $`z`$, the reason for which is quite apparent: The angular diameter $`\theta _{\mathrm{vir}}`$ and the integrated Comptonisation $`\mathrm{{\rm Y}}`$ decreases because of the increasing angular diameter distance $`d_\mathrm{A}(z)`$ that enters $`\theta _{\mathrm{vir}}`$ linearly and $`\mathrm{{\rm Y}}`$ quadratically. Furthermore, clusters accrete matter during their formation history and thus are on average more massive at later times, i.e. at smaller redshifts $`z`$ (see, e.g., Wechsler et al., 2002; van den Bosch, 2002; Zhao et al., 2003). From the physical point of view, the dependence of $`a`$ and $`c`$ on redshift $`z`$ is far from trivial, and therefore, their functional behaviour is described by an empirical approach. Among others, the exponential function provides a good fit to the data, as illustrated by Figs. 10 and 11:

$$x(z)=x_1\mathrm{exp}\left(\frac{z}{x_2}\right)+x_3\text{, where }x\{a,c,s\}\text{.}$$
(37)
The optimised parameters $`x_i`$, $`i\{1,2,3\}`$, for $`x\{a,c,s\}`$ in the gauge functions eqn. (37) are given in table 2 for the case $`q=3`$. It should be emphasised that the parameters stated are only valid for image analysis with the sym2-wavelet, where the maps have been smoothed with a Gaussian kernel with $`1\mathrm{}`$ (FWHM) and the considered cluster sample, which is defined by the selection criteria laid down in Sect. 4.3 and the minimal mass of $`510^{13}M_{\mathrm{}}/h`$.

### 5.8 Influence of CMB noise

Clusters at high redshift $`z`$ are characterised by their small angular scale on which the underlying CMB is represented by a smooth gradient due to Silk damping (Silk, 1968). In this case the wavelet analysis produces the same results irrespective of the CMB noise owing to the distinct morphological feature of the cluster on top of the smooth CMB gradient. Once clusters at lower redshifts reach angular sizes comparable to characteristic scales of CMB fluctuations, the wavelet analysis has to be made more sophisticated. This complication in the wavelet analysis arises because wavelets are primarily suited for determining morphological features rather than solely singling out high amplitude characteristics. Because the angular scale of the clusters ranges between $`10\mathrm{}`$ and $`1\mathrm{}`$, which corresponds to multipole orders of $`\mathrm{}10^3\mathrm{}10^4`$, it suffices to consider the Silk damping tail of the angular power spectrum of the CMB. In the wavelet spectrum $`X_q(\sigma )`$ this translates into an additional approximate power-law component $`X_q^{\mathrm{CMB}}(\sigma )`$, as can be seen from Fig.12:

$$\mathrm{ln}X_q^{\mathrm{CMB}}(\sigma )a_{\mathrm{CMB}}+s_{\mathrm{CMB}}\mathrm{ln}\sigma .$$
(38)
This is due to the discrete sampling of the wavelet moments as well as the inherent statistics of the wavelet spectra of order $`q`$ which can be interpreted as suitably weighted $`q`$-point correlation functions in Fourier analysis (compare Sect. 3.3).

Fig. 13 shows the probability distribution function $`p(s_{\mathrm{CMB}})\mathrm{d}s_{\mathrm{CMB}}`$ of the slopes $`s_{\mathrm{CMB}}`$ following from linear fits to the wavelet moments $`X_q(\sigma )`$ for the range of $`q^{}`$s considered here. Again, the sym2-wavelet was chosen as the analysing wavelet. The slopes $`s_{\mathrm{CMB}}`$ are not well confined, keeping the vast range of angular scales in mind, which in turn will make it difficult to subtract the CMB-contribution to the wavelet spectrum of the combined map.

In order to disentangle the contributions from the CMB noise from those of the cluster, one may pursue different approaches: Among others, CMB fluctuations underneath the cluster can be reconstructed with spline polynomials and successively subtracted. Here, we have masked the cluster and fitted $`5^{\mathrm{th}}`$-order polynomials to the remaining data points. Because the $`y`$-maps and the realisations of the CMB are to leading order combined linearly and because the CMB is a smoothly varying field, it is possible to reconstruct the CMB fluctuations from the environment of the cluster and interpolate to the cluster centre. The reconstructed CMB field can be subtracted from the inital image and by applying wavelet decomposition to the cleaned field one obtains a wavelet spectrum, from which the parameters $`a`$, $`c`$ and $`s`$ can be reliably derived.

### 5.9 Redshift estimation

In order to assess the accuracy of the redshift measurement, a maximum likelihood estimation is performed. The likelihood function is defined as:

$$(z)=\frac{1}{(2\pi )^{3/2}\sigma _a\sigma _c\sigma _s}\mathrm{exp}\left(\underset{x\{a,c,s\}}{}\frac{1}{N}\underset{i=1}{\overset{N}{}}\frac{(x_ix(z))^2}{2\sigma _x^2}\right)\text{,}$$
(39)
and was evaluated for each bin separately, i.e. the index $`i`$ enumerates clusters within the redshift bin under consideration. $`N=30`$ denotes the number of clusters within a single redshift bin. From the position of the maximum in $`(z)`$ the most probable redshift estimate $`z`$ was derived and the accuracy of the estimate followed from the percentiles corresponding to $`1\text{-}\sigma `$ confidence intervals. Fig. 14 shows the estimated redshift versus the real redshift for the cluster sample derived by using all of the three parameters $`a`$, $`c`$ and $`s`$. In comparision, the error bars have become larger by a factor of $`1.5`$ when including the fluctuating CMB, as illustrated by Fig. 15. The measurement is unbiased and the error relative to $`1+z`$ rises slightly with increasing redshift $`z`$.

The results for different analysing wavelets as a function of wavelet moment order $`q`$ are summarised in tables 3 and 4. Clearly, the method starts to fail at redshifts exceeding $`\stackrel{>}{}1`$, when the angular diameter distance $`d_\mathrm{A}(z)`$ develops a plateau and does not cause clusters to appear smaller. The average attainable accuracy is stated relative to $`1+z`$ in order to facilitate comparison to photometric redshifts. The accuracy slightly degrades with increasing $`q`$, which is due to suppression of small wavelet expansion coefficients especially at small scales and the resulting inaccuracy of the fitting formula eqn. (33) used to extract the spectral parameters $`a`$, $`c`$ and $`s`$ from the wavelet spectrum.

Inclusion of the CMB in order to test the applicability of determining morphological redshifts in the case of single-frequency interferometers results in a deterioration of the redshift estimation accuracy of a factor close to 1.5, which is caused by imperfections of the CMB removal by $`5^{th}`$-order spline polynomials.

## 6 Systematics

SZ clusters would be self-similar and would perfectly follow scaling relations provided several requirements are fulfilled: (i) virial equilibrium ($`TM^{2/3}`$), (ii) structural identity, expressed in equal form factors, (iii) a universal baryon fraction and (iv) the absence of heating and cooling processes. Each of these assumptions may be challenged and leads to deviations from the self-similar scaling relations. While the first two points are included in the numerical simulation and are limited by the selection criteria, they increase the scatter in the relations between virial quantities, or equally, the wavelet parameters $`a`$, $`c`$ and $`s`$. Systematic trends of the baryon fraction with cluster mass (see Sect. 6.1) and the formation of cooling flows (Sect. 6.2) need to be assessed separately.

### 6.1 Influence of tilted scaling relations

Analyses of X-ray observations carried out by Arnaud & Evrard (1999) and Mohr et al. (1999) suggest a weak trend of the clusters baryon fraction with cluster mass $`M`$ and a deviation from the universal value $`f_B=\mathrm{\Omega }_b/\mathrm{\Omega }_m`$, which is due to feedback processes like galactic winds that more effectively deplete the ICM of baryons in low temperature clusters compared to high temperature clusters. The variation of the gas fraction, however, is highly model dependent and reflects the mass and density estimation technique: Estimating masses by means of the virial theorem at fixed density contrast, yields a scaling of $`TM^{2/3}`$ and requires the baryon fraction to increase with temperature. Alternatively, using the fit of a $`\beta `$-profile to the cluster, one finds the baryon fraction $`f_B`$ to be nearly constant while a steeping of the slope of the $`M`$-$`T`$-relation is found.

In either of the above mentioned cases the dependence especially of the wavelet parameter $`a`$, which is a logarithmic measure of the SZ flux $`\mathrm{{\rm Y}}`$ would be increased in more massive clusters and would thus increase the scatter in $`a`$ of a cluster sample at a given redshift. The quoted analyses of X-ray data find the baryon fraction to show a relative variation amounting to $``$ 10% at fixed temperature, i.e. at fixed depth of the potential well for a sample of local clusters. Apartfrom the the systematic component, that can in principle beremoved, once high quality X-ray data will improve our understanding of this phenomenon and allows proper modelling, the stochastic contribution can only be constrained to be at most of equal relative influence to $`\mathrm{\Delta }\mathrm{{\rm Y}}/\mathrm{{\rm Y}}`$ as the scatter in morphology.

### 6.2 Cooling flow clusters

In order to estimate the accuracy of the method outlined above, so far we only used adiabatic hydrodynamical simulations which lack of cooling processes. Thus we need to address the influence of cooling flow clusters on our proposed method. After an analytical investigation following Sect. 3.2.1 we compare clusters with and without cooling flows and show how the morphological changes in cooling flow clusters impacts on the wavelet spectra.

#### 6.2.1 Analytical wavelet transform of cooling flow clusters

Instead of a single King profile we assume that the SZ emission of a cooling flow cluster can be described by a double King profile for reasons of analytical feasibility:

$$y(𝒙)=y(r)=\underset{i=1}{\overset{2}{}}y_i\left[1+\left(\frac{r}{r_i}\right)^2\right]^1,$$
(40)
where the second term describes the additional enhancement owing to the cooling flow. Deprojecting this two-dimensional profile by means of Pfrommer & Enßlin (2003) yields:

$$p_\mathrm{e}(R)=\frac{1}{2}n_\mathrm{e}(R)kT_\mathrm{e}(R)=\frac{m_\mathrm{e}c^2}{\sigma _\mathrm{T}}\underset{i=1}{\overset{2}{}}\frac{y_i}{2\pi r_i}\frac{(\frac{1}{2},\frac{3}{2})}{\left(1+R^2/r_i^2\right)^{3/2}},$$
(41)
where $`R`$ denotes the three-dimensional radius and $`(a,b)`$ denotes the $`\beta `$-function (Abramowitz & Stegun, 1965). Thus we obtain for the ratio of the central values of the Comptonisation parameters $`y_i`$

$$\frac{y_2}{y_1}=\frac{p_2r_2}{p_1r_1}\frac{1}{2},$$
(42)
where we inserted typical values for cooling flow clusters, $`p_2/p_13`$ and $`r_2/r_11/6`$. The second order wavelet moment of cooling flow clusters can be obtained in analogy to the non-cooling flow case:

$$X_2^{\mathrm{CF}}(\sigma )=2\pi dkk^5\mathrm{exp}(k^2\sigma ^2)\left|y_1r_1^2K_0(kr_1)+y_2r_2^2K_0(kr_2)\right|^2.$$
(43)
This second order wavelet moment shows an increasing amplitude and a decreasing cutoff parameter compared to the one without a cooling flow.

#### 6.2.2 Numerical analysis

In order to scrutinise these findings we apply our method to adiabatically simulated clusters to which we add an enhanced emission to mimic the SZ emission of the cooling flow. In Fig. 16 the resulting spectra of wavelet moments are shown together with the fitting formula eqn. (33) for increasing wavelet moment order $`q`$.

It can clearly be seen in Fig. 16 that the enhanced emission due to the cooling flow yields a slightly higher amplitude of the wavelet spectrum on small scales. Extracting information from the wavelet spectrum by means of eqn. (33) reveals slightly higher values for the amplitude $`a`$ and smaller values for the cutoff $`c`$ on the percent level. However, this influence is minimised when considering finite instrumental resolution particularly for high redshift clusters. In any case, if a prominent cooling flow is sufficiently well resolved it could be masked and replaced by an interpolation in between the mask boundaries.

## 7 Redshift estimation in a nutshell

This section shall provide a short summary of how to apply our method to an SZ survey for estimating redshifts provided a temperature map of a patch on the sky with resolved images of clusters.

* Once a cluster candidate has been localised at a particular position of the map this cluster and its ambient field has to be cut out. If the number of grid points amounts below $`64^2`$ sampling points, the mesh should be refined by interpolation in order to reach dynamical range of approximately two decades. This is important in order to provide a sufficiently broad range of scales to be probed by the wavelet decomposition.
* The wavelet spectrum of the map is obtained by wavelet transforming the map preferentially using the symlet basis functions (compare Sect. 3.1). The morphological information contained within the wavelet spectrum can be extracted by means of the model function of eqn. (33). In the case of single-frequency observations the ambient CMB field cannot be separated from the SZ signal of a cluster. The method described in Sect. 5.8 might be applied in order to reconstruct the wavelet spectrum of the pure SZ cluster signal.
* The redshift dependence of the wavelet parameters (amplitude $`a`$, cutoff $`c`$, and slope $`s`$) follows the functional form of eqn. (37). However, the single model parameters depend on the definitions of the particular wavelets and the details of the survey, including different sources of noise and the cluster detection criteria. The most promising way of determining the parameters of the gauge functions laid down in eqn. (37) would be to derive them from a learning set of clusters with known (photometric) redshifts. The final redshift estimate of the cluster is most conveniently determined by means of maximum likelihood analysis, as described by eqn. (39).

## 8 Summary

In this paper, a method of estimating the redshift of a cluster based on the wavelet decomposition of its resolved SZ morphology is presented. From a fit to the spectrum of wavelet moments three spectral parameters are derived, that in turn are non-degenerate and indicative of cluster distance. These parameters are utilized, through a maximum likelihood technique, for estimating the cluster’s redshift. In the maximum likelihood technique, empirical gauge functions describing the wavelet parameter’s $`z`$-dependence are used.

First, the method was tested on a simple analytical case: The spectrum of Mexican-hat wavelet moments can be derived analytically for a King-profile, which is known to describe the Compton-$`y`$ amplitude of clusters well. The spectrum of wavelet coefficients as a function of wavelet scale $`\sigma `$, exhibits a break at the cluster scale $`r_c`$ and may thus serve as a measure of the cluster’s size. Additionally, the asymptotic behaviour of the wavelet spectrum in the limit of $`\sigma r_c`$ and $`\sigma r_c`$ can be understood. The derivation of wavelet moments of order $`q=2`$ is analogous to considering the Fourier power spectrum of the Compton-$`y`$ map, filtered with Fourier transformed wavelet. The shape of the spectrum of wavelet moments of order $`q=2`$ from the analytic calculation is consistent with one obtained by applying wavelet decomposition to simulated SZ cluster maps.

The method was then applied to set of numerically simulated SZ clusters with $`1\mathrm{}`$ (FWHM) resolution – comparable to the resolution of future SZ experiments. The sample comprises 690 cluster maps distributed in 23 redshift bins, which is a comparably large cluster sample. The clusters are chosen such that they are not in a merging state and their SZ image is not too elongated, two criteria that favour clusters close to virialisation. Additionally, in order to simulate single-frequency observations, the cluster maps were combined with realisations of the CMB that constitute the main source of noise.

The method was tested for a range of wavelet functions (e.g., symlet, coiflet, Daubechies, biorthogonal). The average attainable accuracy in estimating redshifts is found to be almost independent of the specific functional form used, although the symlet basis yielded the best results. However, the method could benefit from improvements concerning the choice of the wavelet basis. For instance, one could try to construct an optimised wavelet specifically for $`\beta `$-profiles, that yields maximised wavelet coefficients $`\chi (𝝁,\sigma )`$.

As expected, there is only a weak change in accuracy with respect to the order $`q`$ of the chosen wavelet moment $`X_q(\sigma )`$. This, however, is most likely to change when applying the wavelet analysis to noisy images, because for increasing choices of $`q`$, uncorrelated noise is suppressed relative to the cluster’s signal and concentrating on higher values for $`q`$ should provide a more robust measurement of the set of structural parameters $`a`$, $`c`$ and $`s`$. The increment of $`q`$ itself is limited by numerics – this is the case when the wavelet moment $`X_q(\sigma )`$ is dominated by the largest wavelet expansion coefficient $`\chi (𝝁,\sigma )`$, and does not reflect anymore the dependence on the wavelet scale $`\sigma `$. In this limit, the wavelet spectrum would exhibit a generic power law behaviour: $`X_q(\sigma )\sigma ^{\gamma (q)}`$ for large $`q`$. The structural parameters $`a`$, $`c`$ and $`s`$ were found to depend on redshift $`z`$ by a simple exponential (eqn. (37)). The free parameters in this equation can be determined from a (relatively small) sample of SZ cluster images with known redshift.

The accuracy of determining cluster distances has been assessed by maximum likelihood estimation. The method yields accuracies of $`45\%`$ relative to $`1+z`$, which is competitive with photometric redshifts, but reaches out to larger distances. At redshifts exceeding $`z\stackrel{>}{}1`$, the accuracy is expected to degrade because the angular diameter distance $`d_\mathrm{A}(z)`$ starts to level off and thus sets the limit of applicability. For single frequency data, the CMB fluctuations can be removed with a simple polynomial reconstruction approach; the accuracy in the redshift estimation is then decreased to $`67\%`$.

In this work we have considered two major systematic effects that might degrade the accuracy of the method. The first is the varying baryon fraction with cluster mass, which has been studied only for local cluster samples. While the systematic trend could in principle be corrected for, the stochastic contribution will always add to the uncertainty of the distance determination. Another systematic is the influence of cooling flows at the cluster’s centre. In this case we have been able to show that the uncertainty it adds to the redshift estimate is very small, mainly because the volume occupied by the cooling flow region is limited to the cluster’s core.

Although the result in the distance estimation is stated in terms of redshift, it should be emphasised that a specific cosmology is assumed, which is needed for converting the observables, namely the wavelet parameters, to a distance estimate. The distances following from the analysis have been expressed as redshifts because of their elementary interpretation, but the implicit assumption of an underlying cosmology should be kept in mind when comparing to e.g. photometric redshifts. For that reason, the precision of the method presented is limited by the accuracy to which the cosmological parameters are known. Apart from being a distance indicator, the redshift also plays the role of an evolutionary parameter.

Comparing this work to the pioneering paper by Diego et al. (2003), our expectations concerning the accuracy of morphological redshifts are even more optimistic: Without fitting $`\beta `$-profiles to the observational data, it is possible to describe the cluster’s SZ morphology by solely relying on wavelet decomposition. Also, we describe the spectrum of wavelet moments with a small set of structural parameters, that have a lucid physical interpretation, provide a non-degenerate distance measurement and enable redshift determination owing to their monotonic decline with redshift. Furthermore, the redshift dependence of the structural parameters is calibrated with the data set itself without relying on prior and simplifying assumptions. In spite of the small number of observables considered here, the accuracy in the redshift estimation of this method is doubled, in comparison with Diego et al. (2003), even for a single frequency experiment.

## Acknowledgements

The authors would like to thank Volker Springel and Lars Hernquist who kindly provided us with their numerical simulations and Antonaldo Diaferio, Matthias Bartelmann and Simon D.M. White for clarifying discussions and many useful comments.


format: text

text:
# In the search for the low-complexity sequences in prokaryotic and eukaryotic genomes: how to derive a coherent picture from global and local entropy measures

## I Introduction

In the recent past there has been a significant interest in the search for both long- and short-range correlation in DNA sequences russi ; voss ; stanley ), The main results of these papers have been that the amount of long-range correlation increases when moving from prokaryotes to eukaryotes and from coding to non-coding sequences. This paper is devoted to discussing the same issue using two new theoretical tools that were not yet available to the authors of the earlier papers. The former method, called Diffusion Entropy (DE), was developed by the authors of Refs.giacomo ; nicola ; PRE  for the purpose of defining the asymptotic properties generated by the long-range correlation. When this method is applied to DNA sequences, it affords global information, this being unambiguous only in the ideal stationary case. The DNA sequences are non stationary and the adoption of the DE method must be supplemented by a non-stationarity indicator. This important indicator has been proposed in Ref. giulia1 . We refer to it as Non-Stationarity Entropic Index (NSEI) method. This is so because, as we shall see, this method is entropic in nature, it vanishes in the stationary case, and it is larger the larger the strength of non-stationarity. The measurement of long-range correlation is a delicate issue. The detection of non-stationarity is another delicate issue. The NSEI method addresses properly an even more delicate issue, this being the non-stationary character of long-range correlation.

With the use of these two brand new techniques we prove that moving from coding to non-coding and from eukaryotes to prokaryotes both the non-stationarity degree and the amount of long-range correlation increase. The deviations from stationarity and randomness are a sign of the role played by selection during the process of life evolution. The same argument applies to the increasing amount of non-stationarity, since this phenomenon is thought to be associated with the selection-induced increasing variety of molecular functions.

To make more solid this important biological conclusion we proceed as follows. We take into proper account the ideas that have been developed in the last few years to properly model the origin and the functions of the constraints on the DNA sequences randomness ruffodna . More specifically, we shed light into this issue by analyzing in the light of life evolution some key Low-Complexity Sequences (LCS). These LCS are defined as sequences containing only AT (homo-weak) or GC (homo-strong), or purines (GA), or pyrimidines (CT). The reason for choosing these particular non-random sequences lies in the fact that their presence is known to affect the DNA local conformational landscape, giving rise to non-B structures (the B structure is the usual double helix), and/or affecting the degree of bending, curving, twisting and rolling of DNA molecules. Conformational features are known to be true 3-dimensional codes for the interaction of DNA with proteins, a necessary condition for its involvement in chromosome organization and in the key functions of transcription and replication buiattiinpress .

The outline of the paper is as follows: In section I.1 we describe some basic properties of LCS distributions and their influence on DNA structure. In section II we describe the global analysis trough the DE of different genomes and compare it with LCS distribution data, while in section III we illustrate the NSEI method. Section IV aims at showing the results of two methods at work, namely, the NSEI and a neural network trained by means of LCS lengths to distinguish between coding and non-coding sequences. Finally Section V is devoted to drawing some conclusions.

### I.1 Distribution of LCS and their putative functional meaning

Provata and Almirantis almirantis  analyzed the size distribution of purine and pyrimidine clusters and found that coding and non-coding sequences yield an exponential decay and a power law decay, respectively. We also note that in eukaryotes the coding sequences relative amount is very low, whereas it is the largely dominant fraction of prokaryotes genomes. Both properties suggest to make an analysis of LCS length distributions in species placed in different positions in the “tree of life”buiattiinpress . In a preliminary paper buiatti2002  an algorithm was used to evaluate the LCS density in a large number of bacterial, archaeal and eukaryotic genomes. A clear difference was observed in the distributions of LCS relative amount between coding and non-coding regions and between prokaryotes and eukaryotes. In all genomes non-coding and Eukaryote sequences show the highest amount of LCS. To understand whether high-density values were accompanied by a random dispersion or by the presence of long LCS stretches, we decided to carry out an analysis of LCS length distributions. Herein the attention is centered on LCS known to give rise to non-B DNA conformations and to be involved in the control of gene expression and recombination, namely the homo-weak and the homo-strong. In the case of AT-rich (homo-weak) sequences we found higher values of LCS lengths in non-coding sequences vs. coding sequences and in eukariotes vs. Prokaryotes.

Fig. 1 shows this property for homo-weak noncoding sequences in some genomes, representative of archaea (M. Jannaschii), prokaryotes (B. Subtilis), unicellular eukaryotes (S. Cerevisiae) and multicellular eukaryotes (H. Sapiens, Chromosome I). We see an increasingly larger deviation from a single exponential distribution going from prokaryotes to more complex organisms with a larger amount of non-coding sequences. This property is less trivial than it may appear, since it is true also for GC-rich (homo-strong) sequences, as denoted by the curves with white symbols in Fig. 1. All these curves show a clear deviation from Poissonian behavior, even if their inverse-power-law decay has a large negative exponent (about $`5`$). However, for all genomes, the curves for homo-weak LCS have a larger cut-off value than the corresponding ones for homo-strong LCS. This fact is indicative of a selection pressure against the lengthening of homo-strong sequence possibly to be attributed to their rigidity and conformation landscapes different from those of long AT-rich tracts. This leads us to choose homo-weak occurrence as a marker in the correlation analysis, as it will be explained below. Indeed, the presence of a decay slower than Poissonian suggests the possibility that correlation with an extend long range can be triggered by the presence of homo-weak and homo-strong LCS.

Remarkably, the authors of pomodori  performed a thorough study of intra and inter specific variability in upstream and coding sequences of a series of key genes in Lycopersicum (tomato) (the gene for ACC synthase), Nicotiana (the gene for Phytochrome), Primate species (the genes EMX and OTX 2 involved in brain development). All data showed a far higher variability in non-coding than in coding sequences, preferentially localized in LCS. In other words, different varieties have different lengths of AT-rich segments in the non-coding region before the gene. Quantitative RT-PCR data (Reverse Transcription-Polymerase Chain Reaction), a method to establish the amount of messenger-RNA, showed a clear variation of gene expression corresponding to different lengths of the homo-weak sequences.

Remarkably, the authors of pomodori  performed a thorough study of intra and inter-specific variability in ACC-synthase upstream sequences. They studied ACC-synthase upstream sequence in Nicotiana spp phytochrome regulatory sequences and in the same regions for the primates genes EMX-2, OTX-2. The results clearly showed that these LCS tend to be hyper-variables, i.e. they show a mutation frequency significantly higher than other DNA stretches, moving from and individual to another of the same species as well as from species to species, coherently with the aforementioned data on recombination. The interesting conclusion is that the three-dimensional code regulating the gene expression may be associated to a strongly dynamical DNA structure, with statistical constraints. This structure is probably associated to a process with long-range correlation. In the following section we shall analyze the long-range correlation present in the genome sequences and the relation between this property and the presence of LCS.

## II The diffusion entropy method applied to homogeneous sequences

### II.1 Global Diffusion Entropy

In recent papers giacomo ; nicola ; PRE  a new way of revealing long-range correlation, based on the detection of anomalous scaling has been proposed. In short, one defines a “marker” on a time sequences, and evaluates the probability $`p(x;t)`$ of having the number $`x`$ of markers in a window of length $`t`$. The evaluation of $`p(x;t)`$ is done by moving a window of length $`t`$ along the sequences and counting how many times one finds $`x`$ markers inside this window. At this stage $`p(x;t)`$ is obtained by dividing this number by the total number of windows of size $`t`$, which is obviously $`Nt+1`$, where $`N`$ is the total length of the sequence. In the case of large values of $`x`$ and $`t`$, the continuous approximation is legitimate, and, in the ergodic and stationary condition, a scaling relation is expected, namely In case of large values of $`x`$ and $`t`$, the continuous approximation is legitimate, and, in the ergodic and stationary condition, we expect scaling, namely, we expect

$$p(x;t)=\frac{1}{t^\delta }F\left(\frac{x}{t^\delta }\right),$$
(1)
where $`\delta `$ is said to be the scaling index and $`F`$ is a function, sometimes called “master curve”. If $`F`$ is the Gauss function, $`\delta `$ is the known Hurst index, and if the further condition $`\delta =0.5`$ is obeyed, then the process is said to be Poissonian, namely there is no long-range memory regulating the occurrence of markers in time. It is straightforward to show that the Shannon Information

$$S(t)=_{\mathrm{}}^{\mathrm{}}𝑑xp(x;t)\mathrm{ln}p(x;t)$$
(2)
with conditiom (1) leads to

$$S(t)=k+\delta \mathrm{ln}t,$$
(3)
where $`k`$ is a constant. A linear fit of $`S(t)`$ in log-normal paper allows an extimation of $`\delta `$.

### II.2 The CMM model of DNA

In recent years many groups have agreed to consider DNA sequences as a mixture of long- and short- range correlation, the latter blurring the strength of the former. A simple model, originally proposed by Araujo et al. araujo , and later adopted independently by the authors of Refs. maria , marcobuiatti2  and zanettedna , to name a few, is the Copying Mistakes Map (CMM)maria . This is the superposition of two models, one consisting of a pure random choice, and the other of an intermittent generator of homogenous sequences, yielding long-range correlation. Both models provide a sequence, where the $`i`$-t nucleotide of the resulting DNA sequence is taken from the first component (the white one) with probability $`1ϵ`$ and from the second with probability $`ϵ`$. The resulting correlation function is proven maria  to be

$$C(t)=\delta _{i,0}+ϵ^2\left(\frac{T}{T+t}\right)^\beta ,$$
(4)
where $`\delta _{i,j}`$ is the Kroeneker delta, while $`T`$ and $`\beta `$ are two positive parameters of the intermittent model. Long-range memory is normally characterized by the condition $`\beta <1`$, which makes the correlation function (4) not integrable.

As shown in ref. maria , the second moment and the Hurst analysis reveal a short-time behavior which is dominated by the random component, while the correlated part dominates in the long-time limit. By the same token, DE method applied to CMM sequences, yields for $`S(t)`$ a curve that starts with a slope of $`\delta =0.5`$, and then after a knee, tends asymptotically to higher value of $`\delta `$ corresponding to the Lévy scaling $`\delta =1/(\beta +1)`$. This behavior, as reported in Ref. latora , is illustrated in fig.2, which shows that the position of the knee is a monotonic function with respect to $`ϵ`$: the larger $`ϵ`$ the smaller the position of the knee.

### II.3 DE at work on real genomes

As earlier shown, the DE method establishes the values of two quantities: the ratio of the correlated to the uncorrelated component intensity of the CMM model (through the position of the knee) and the anomalous scaling index $`\delta `$. As an example, we show in fig. 3 the DE analysis on the chromosome IV of Saccharomices Cerevisiae. In Tab.1 we report the results corresponding to different genomes grouped according to whether they are bacteria, archaea or eukaryotes. All these analyses were performed using as marker the homo-weak occurrence in the series. This means that we put $`1`$ at position $`i`$ if the $`i`$-th nucleotide and the previous one are weak bases (AT), and $`0`$ elsewhere. This choice was suggested by the non-exponential decay in the length distribution of homo-weak sequences shown in fig.1a. We made the hypothesis that this non-exponential behavior can be connected with the long-range correlation observed in the genomes.

Tab.1 shows two interesting results: the former is that the value of $`\delta `$ does not seem to depend on the position in the evolutionary tree. The latter is that the position of the knee seems to change with moving from bacteria to archaea and eukaryotes. As shown in the table, the weight of the correlated component in bacteria is larger than in eukaryotes, archaea being somewhat intermediate between bacteria and eukaryotes. It is also worth noticing that within eukaryotes, the unicellular, almost intronless, S. Cerevisiae has the largest knee value.

## III The Non-Stationarity Entropic Index

In order to study the local properties of DNA sequences with special focus on the distribution and variability of homogeneous sequences, we use a new method, the NSEI method. This method is derived from the DE method giulia1  as an earlier method, called CASSANDRA algorithm giulia1 ; granada . The DE method affords global information, while both CASSANDRA and NSEI method aim at affording local information. The CASSANDRA method measures the rate of transition from dynamics to thermodynamics through comparison with ordinary Brownian motion. The NSEI method focuses on the local deviation from the stationary behavior. The first step of the NSEI method is to build the diffusion entropy function $`S(t)`$ for each window of size $`L`$ with $`Lt`$, thought of a complete sequence. In this way we obtain several functions $`S_j(t)`$, where the subscript $`j`$ refers to the portion of the sequence used to calculate $`S_j(t)`$. Formally, we write

$$x_j(t,l)\underset{i=j+l}{\overset{j+l+t}{}}\xi _i,\mathrm{with}\mathrm{\hspace{0.25em}\hspace{0.25em}0}<l<Lt$$
(5)
where $`\xi _i`$ is the DNA sequences after marking the homogeneous sequences, as earlier ilustrated, and $`j`$ and $`L`$ denote the beginning and the length of the big window, respectively. We assign to the integers $`j`$ multiple values of a given integer $`J`$, namely $`j=0,J,2J,\mathrm{}`$, thus avoiding too large overlaps between two windows of length L, corresponding to two consecutive choices of $`j`$. At fixed $`t`$ and $`L`$ we obtain for different $`l`$ a probability distribution $`p_j(x,t)`$. Calculating the entropy of this distribution, we obtain

$$S_j(t)dxp_j(x,t)\mathrm{log}p_j(x,t)$$
(6)
Now we are ready to define the non- stationarity index, indicated as $`\mathrm{\Xi }_j`$, as follows

$$\mathrm{\Xi }_j=\underset{0}{\overset{t}{}}dt^{}\left(S_j(t^{})\overline{S}_j(t^{})\right)$$
(7)
where

$$\overline{S}_j(t)\frac{1}{m}\underset{i=jm1}{\overset{j1}{}}S_i(t)$$
(8)
Note that $`m`$ is the number of windows involved by the average process of Eq.(8). Depending on the step $`J`$ by which $`j`$ increases we obtain that the comparison between local and previous diffusion behavior involves $`mJ`$ points of the sequence. Note that if the sequence is perfectly stationary we have $`\mathrm{\Xi }_j=0`$, otherwise we have a value different from zero.

## IV DNA function and local correlation

In Section I.1 we saw that non-coding sequences contain a number of homogeneous sequences larger than coding sequences. Experimental evidence shows that LCS tend to occur either in the upstream region of each gene, i.e. immediately before the occurrence of a gene, or in introns. It is also known that both upstream sequences and, to a lesser extent, introns, are the regions endowed with regulatory roles in which LCS are involved. In this section we evaluate the change of local correlation, according to the NSEI method, and we see that this change is often associated with the beginning of a coding portion. Thus, using biological consideration, we are led to establish a plausible link between the change of DNA sequence correlation and the functional change of the sequence. In other words, we expect that the boundary between a coding and a non-coding segment is revealed by either a positive or negative peak of the NSEI index, as a function of the position $`j`$ along the sequence. In Fig. 4 we show an example where this hypothesis is verified to a satisfactory extent. As earlier stated, the analysis was carried out using as marker the occurrence of homo-weak sequences. Peak recognition has been realized in the following way. The data have been smoothed through a moving linear interpolation with ten points. The maximum indexes were identified looking at the “derivative” thus obtained; only maximum indices with a value larger than the standard deviation, calculated over the whole signal, were selected.

The NSEI method should reveal the beginning and the end of a LCS with peaks. It is important to point out that the NSEI method, although being the most convenient indicator of non-stationarity, within the theoretical framework of DE method, does not have a perfect efficiency. This means that the NSEI signals correctly the existence of a boundary between two distinct regions. This is a true positive, denoted by the symbol $`TP`$. There are cases, however, when the NSEI signals a boundary that does not exist. This is called a false positive, and it is denoted by the symbol $`FP`$. By the same token we denote by $`TN`$ and $`FN`$, the true and false case, respectively, of a missing boundary. The total number of borders is given by $`FN+TP`$, while the total number of cases signaled by the NSEI, true or false, is given by $`TP+FP`$. It is worth noticing that the observed peaks, in addition to corresponding to true peaks, must be found in position close enough to the position of the real borders, namely at a distance smaller than a critical value $`T`$, which is of the same order as $`L`$.

We are now in the proper position to establish the accuracy of the NSEI. We can judge this method satisfactory if we prove that it affords a number of correct guesses significantly larger than the random choice. To assess this issue we apply to the NSEI method a technique of assessment frequently used in Information Extraction precall . We define two quantities, Precision and Recall. Precision $`P`$ is a measure of how accurately the displacements of coding/non-coding boundaries are guessed, while Recall $`R`$ is a measure of how many of these regions are correctly located.Formally, we write

$$P=\frac{TP}{TP+FP}$$
(9)
and

$$R=\frac{TP}{TP+FN}$$
(10)
We compare the $`P`$ and $`R`$ of our predictor to a baseline score associated with $`TP+FP`$ random choices for the position of the boundaries. The “true” classification was done according to documents accompanying DNA sequences in Genbank database, regardless whether the annotation was experimental or putative. Although this makes our esteem or $`P`$ and $`R`$ questionable, we think that the improvement with respect to the random baseline will not diminished after a more accurate experimental annotation. Results are summarized in table 2. As shown in the table, the NSEI method reveals statistically significant improvement with respect to the baseline, especially for eukaryotes, where the marker adopted (homo-weak) has a larger biological significance.

We end this section by reporting results stemming from a completely different way of detecting coding and non-coding sequences. This alternative way confirms the correlation between the changes of the amount of LCS and entropy changes, when moving from the coding to the non coding regions. A three layer feed forward neural network was built and trained on homogeneous sequences (homo AT; homo GC; homopurine; homopyrimidine) neuralnets . The network involves eight neurons in the input layer, four in the intermediate, one in the output one. The network has been trained on 100 tracks each of 300 nucleotides of known coding and on 100 for known non coding regions. Inputs of the network are 1) percentage of the sequence covered by homogeneous tracks 2) average length of homogeneous tracks in the observed sequence. Then 900 tracks have been screened with the trained network for each class. This operation has been performed for each species listed in the table 3. The precision level of annotation (identification of the two classes) is very high as seen from Q and W values. For true positives (Q), moreover, eukaryotes seem to perform better than prokaryotes, coherently with the higher non-coding/coding ratio.

It is worth noticing that this second way of detecting the change of the DNA function is complementary to the first one. Indeed the NSEI method, which does not require any form of training, yields a model-oriented location of the function change on a scale comparable with $`L`$. On the other hand the neural network technique is more precise and can become very efficient when performed on the basis of evidence obtained by means with NSEI.

## V Conclusions

In this paper we discussed how to derive a coherent picture from the study of LCS, local properties, and long-range correlation, global properties. Thanks to the DE technique we showed that long-range correlation seems to have a quite homogeneous behavior when looked at using a large scale. Variations in the transient behavior, however, confirm an increase of the non-random component with moving from bacteria to archaea, from archaea to the unicellular eukaryote S. cerevisiae, and from this to multi-cellular plants and animals. This increase throughout the evolution scale is consistent with that of non-coding sequences, which are known to cover less than 10% of prokaryotes genomes and well over 90% of the eukaryotes ones. A plausible hypothesis, to explain these properties in terms of function, is that non-coding sequences can be more freely “filled” with “hidden codes”. These hidden codes are mainly used for the recognition of protein molecules, necessary for the organization and function of the DNA sequence and other ligands. Remarkably, these codes, and in general non-coding sequences as well, tend to be hyper-variable, due to the presence of LCS, suggesting a positive role of mutational noise in regulatory sequences (see Ref buiattiinpress  for a through discussion).

To afford some more details on why to connect LCS, conformational landscapes, regulation and structure of genomes, we should note that some constraints may play a role in eukaryotes for DNA “packaging” in nucleosomes and chromosomes. Moreover, it is well known that DNA activation for transcription and replication is due to the recognition by proteins of DNA regions, and vice-versa, based on complementary conformations leading to formation of complexes. Protein-DNA complexes involve at the same time a high number of molecules where interaction needs a specific global organization favored by DNA curving, bending, twisting and rolling at the right points. This may require the presence of non-random tracts with the needed conformational landscape. The strategy reported herein, namely the adoption of NSEI and of neural networks, fulfills the biologist’s need for research leading to the characterization of the living state of matter (LSM) buiabuia  on one hand, to the identification of specific functional sequences in sequenced genomes on the other.

We investigated on the connection between correlation and LCS, and between LCS and the regulation of a gene transcription. We showed that the variation in the statistics of homo-weak sequences, i.e. a particular kind of LCS, yielding specific DNA conformational landscapes, strongly correlates with changes in DNA function and, especially, that these sequences signal the beginning of a coding region. We studied these variations with the help of the NSEI algorithm obtaining satisfactory results. The accuracy achieved, although affording significant information, is not high, in comparison with other annotation methods currently used by molecular biologists. In fact, with some of these methods, the percentage of failure is now around 10%-15%. Better results were herein obtained with a neural network trained on LCS, the class of non-random sequences. The NSEI method, however, a model-based technique, affords a way to scan the whole genome, in the search for non-stationary segments, where a transition from a given degree of correlation to another occurs: The putative boundaries between coding and non-coding DNA may generate the segmented input tracks for the neural network, a method with a larger accuracy for the final categorization. The NSEI method is based on a DNA dynamic model, compatible with the DE method, which signals the long-range corelation of the sequence. If the correlation strength undergoes local changes, these are signaled by the NSEI method. Moreover, the increase in complexity in the “tree of life” correlates with an increase in non-coding DNA and in the number of coding/non-coding boundaries. A high precision on the same task of recognizing coding sequences has been obtained by other authors with algorithms based on specific motifs known to be present in non-coding regions or exploiting existing periodicities herzelgrosse2000 . Probably an even higher precision could be obtained through the combination of all methods. It is therefore quite important that the NSEI results were obtained without using all these properties but only the variation in the statistics of homo-weak sequences.

A reasonable interpretation of these data all taken together, could be that several kinds of functionally constrained sequences have been fixed in genomes throughout evolution and that average higher correlation values in non-coding sequences are a result of this process. However at a local level, while to some extent LCS can be present also in coding regions, precise annotation (distinction between coding and non-coding sequences) has to rely on the detection of constraints (codes) specific for each class of DNA regions. For instance, while a 3-periodical behavior, in phase with the triplet “universal” code is typical of coding regions maria ; herzelgrosse2000 , specific motifs and rather long LCS have been shown to be distinctive of non-coding regions, which are free from coding constraints maria . In other words, the relatively low level of constraints corresponding to long-range correlations (and to the presence of LCS) in transcribed and translated sequences can be attributed to the fact that in that case selection has not been acting on DNA, but on proteins, and there is not a direct relationship between constraints in the two classes of molecules.

As a final remark, looking at LCS, we point out an almost paradoxical result, namely that hypervariable sequences carry a constraint. LCS hypervariability is known to be favored in regions of genes where variation, leading to protein variability, is essential. For instance, it can be useful for the right response to hypervariable pathogens. Lower complexity levels in eukaryotes seem therefore to derive from a wide series of constraints ranging from periodicities to different kinds of short-range and long-range correlations and to the presence of LCS of varying length. This hypothesis is supported by the present work, suggesting the NSEI as a tool for analyzing data, when looking at this kind of constraints, because of its locality and its wise use of statistics. Furthermore, we want to point out that the NSEI emerges from the theoretical background of Ref. giulia1 . The DE method was originally introduced as an efficient technique to evaluate scaling, this being a global property, implying stationary condition. The NSEI is the best way, known to us in this moment, to address the challenging issue of non-stationarity from the same dynamic approach to complexity as that behind the foundation itself of the DE method. From a conceptual point of view, more than from an application point of view, we want to stress that this technique of analysis, resting on entropy formalism, might have the important role of helping the foundation of the LSM theory buiabuia ; gerardo . As explained in Ref. buiabuia , in fact, the relation between hyper-variable sequences and regulation is a key aspect of the LSM perspective, and the investigation of this issue based on the NSEI might make easier for the advocates of the LSM perspective to express the main ideas of this theoretical proposal with the language of statistical mechanics, probably anomalous statistical mechanics.

*Acknowledgments* PG gratefullly acknowledges financial support from ARO through Grant DAAD19-02-0037