text:
# Xmm–newton observation of Abell 1835 : temperature, mass and gas mass fraction profiles.

## 1 Introduction

Clusters of galaxies are dark matter dominated and thus are ideal objects to obtain information on this kind of material. A straightforward way to determine the mass profile of a cluster, which consists of roughly 80% dark matter is to use the hydrostatic approach based on the hot x–ray emitting intra-cluster medium (hereafter icm). In order to get reliable results with this approach it is necessary that the examined galaxy cluster is in a relaxed state without ongoing merger activity, which might create shock waves and bias the mass estimate. The validity of the hydrostatic approach for relaxed clusters has been tested via hydrodynamic simulations (Schindler 2001 ; Evrard, Metzler & Navarro 1996). To solve the hydrostatic equation one needs to know the density as well as the temperature distribution of the icm. The latter one requires x–ray observatories, which provide at the same time good spatial as well as spectral resolution. X-ray observatories such as Rosat, Asca, and Beppo-Sax only provided limited combinations of these two requirements and thus the determination of cluster temperature profiles based on these data were subject to large error bars and sometimes different studies on the same object with the same data came to very different conclusions (Irwin & Bregman 2000 ; Markevitch et al. 1998). A subject of hot debate was the question whether there exists a universal declining temperature profile in clusters as suggested by Markevitch et al. (1998) or not. Recently launched x–ray satellites namely Chandra (Weisskopf et al. 2000) and Xmm–newton (Jansen et al. 2001) fullfill now all requirements for an accurate determination of the temperature distribution in clusters up to large radii.

We present here a spectro-imaging study and subsequent mass determination analysis of the medium distant (z=0.25) relaxed galaxy cluster Abell 1835 (Allen et al. 1992) observed during the Xmm–newton Performance Verification phase. Abell 1835 is a relaxed cluster, which hosts a cooling flow in its center (Allen et al. 1996). It belongs to the cooling flow clusters for which no multiphase gas down to low temperatures has been found in the Xmm–newton rgs data (Peterson et al. 2001), and which casts serious doubts on current cooling flow models (Fabian et al. 2001, Molendi & Pizzolato 2001).

We determine the total and gas mass profile of the cluster up to 1.7 Mpc which corresponds to 0.75 virial radius (see section 5). The ratio of gas mass over total mass allows us to determine a lower limit for the baryon fraction of the cluster which can be considered as representative of the mean baryon fraction of the Universe (see White et al. 1993 and references therein). Knowing the baryon density of the Universe based on studies of primordial nucleosynthesis we can give constraints on the cosmological matter density parameter $`\mathrm{\Omega }_\mathrm{m}`$.

Our paper is organized as follows : in section 2 we describe the observation. In section 3 we present the data analysis with specific emphasis on vignetting corrections and flare rejection. In section 4, we present first results and a method to correct for different background components : a) particle induced background and b) astrophysical background, which varies across the sky (Snowden et al. 1997). The method is based on a double background subtraction. Section 5 presents our gas and total mass determination and estimate of the baryon fraction of Abell 1835. This is followed in section 6 by the discussion of our results and our conclusions in section 7.

Throughout this article, we use a cosmology with H$`{}_{0}{}^{}=50`$ km/s/Mpc, $`\mathrm{\Omega }_m=1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ (q$`{}_{0}{}^{}=0.5`$). One arc-minute corresponds to 297 kpc at the cluster redshift of 0.25.

## 2 Observation

Xmm–newton observed Abell 1835 in revolution 101 on June 28th 2000. The observation (id 98010101) was taken during the phase of Performance Verification. The total exposure time was 60 ksec. Since we are interested in spatially resolved spectroscopy we concentrate here on data from the epic-instruments (Turner et al. 2001 ; Strüder et al. 2001). During the observation of Abell 1835 the mos1 camera was operating in Large Window mode, which shows a relatively high internal background level, and for which we cannot yet properly account for. Since our work is very sensitive to background uncertainties we discard mos1 data of Abell 1835 in the following and only use data coming from the mos2 and pn camera. For all cameras thin1 filter was used.

## 3 Data processing

We generate calibrated event files using the tasks emchain and epchain in the xmm sas version 5.0 111http://xmm.vilspa.esa.es/sas. For the mos2 data set, we take into account event patterns 0 to 12 and for pn data only single events (pattern=0).

### 3.1 Flare rejection

The light curve of this observation (see figure 1) is not constant and large variations of intensity are visible. These variations, which we will call “flares” in the following, are caused by soft energy protons produced by solar activity, and which are only visible outside the Earth’s magnetosphere, which corresponds to the orbit of Xmm–newton (see Ehle et al. 2001).

Unfortunately, the background induced by flares cannot be easily accounted for since it shows spectral variability. Furthermore, the flares reduce the observed signal-to-noise ratio. Therefore it is best to discard flare periods. This reduces the effective exposure time, but improves the signal-to-noise ratio significantly.

For our analysis we only consider events inside the field-of-view (fov) and ignore all events outside the fov.

Flares show a hard component observable at high energies easily detectable above 10 keV where the mos cameras are hardly sensitive to X-ray photons. We choose the 10 to 12 keV energy band to monitor the background for the mos2 camera. In order to optimize the background rejection, we analyze the relationship between effective exposure time and the threshold we apply. The results are shown in figure 2 for Abell 1835 and Lockman Hole (revolution 70 and 71) observations. We can see that the behaviour of exposure time versus rejection threshold is similar for all observations, with a flattening at a threshold of about 15 counts/100 s. Thus we can say that this is the optimal threshold for flare rejection, since it limits exposure time loss. A higher threshold does not significantly increase the exposure time and a lower threshold limits severely the remaining effective exposure time. Therefore we adopt in the following analysis a threshold of 15 counts/100 s in the 10 to 12 keV energy band. The same procedure is applied for the pn data (see figure 3). The pn camera provides higher quantum efficiency at high energies than the mos detectors. Because of this we shift our energy window for monitoring the background level to the 12 to 14 keV energy band. Our adopted threshold is 22 counts/100 s.

Of course, the adopted threshold criterion for flare rejection must also be applied to the observations chosen to determine the background.

After applying this flare rejection method, the remaining exposure times are 25750 s for mos2 and 25950 s for pn.

### 3.2 Vignetting correction

The collecting area of the x–ray telescopes aboard Xmm–newton is not constant across the field of view. The effective area decreases with increasing off-axis angle and is furthermore photon energy dependent. This effect (vignetting) must be taken into account in the case of extended sources, like galaxy clusters, which cover large fractions of the field of view.

To correct for this vignetting effect, we use the method described by Arnaud et al. (2001a). Briefly, this method consists of calculating a weight parameter for each event. According to event energy, this parameter is defined as the ratio of on–axis effective area to effective area at the event position on the detector. Thus, all weight parameters are greater or equal to unity.

## 4 Data analysis

### 4.1 Morphology

As can be seen in figure 4, Abell 1835 is a compact cluster. There is no strong evidence of substructures which suggests that the cluster is in a relaxed state. The intensity is very peaked in the center due to the presence of a cooling flow (Allen et al. 1996; Schmidt, Allen & Fabian 2001; Peterson et al. 2001).

### 4.2 Spectra

Since we apply the described vignetting correction method, we can use the on–axis response (i.e. effective area $`\times `$ redistribution matrix) files : m2\_thin1v9q20t5r6\_all\_15.rsp for the mos2 camera and epn\_ff20\_sY9\_thin.rsp for the pn. There exist different response matrix files for the pn-chips, which correspond to the different nodes. We use Y9, which is the node farthest away from the readout, with lowest resolution. We tested the different response files of the different nodes and could not find any detectable difference when performing spectral fits. Due to the remaining uncertainties of the response files at very low energy, we exclude events with energy below 0.3 keV for all spectral fits.

#### 4.2.1 Background treatment

Due to the fact that clusters of galaxies are extended and generally low surface brightness sources, the correct determination of the background is important for data fitting and modelling. Since cluster emission decreases with radius from the center, the background component becomes more and more important with increasing radius.

After removing the flare component by excising high background intensity time intervals, two other background components remain :

* high energy particles like cosmic–rays pass through the satellite and deposit a fraction of their energy on the detector. They are not affected by the telescope vignetting. These particles induce instrumental fluorescent lines coming from material aboard the satellite. This kind of background is well monitored by blank–sky observations based on several high galactic latitude observations in which x–ray sources were excised. These observations were compiled by D. Lumb (2002)222They can be retrieved from the Vilspa ftp site : ftp://xmm.vilspa.esa.es/pub/ccf/constituents/extras/background for each epic camera ;
* the astrophysical x–ray background. The soft component of this background is position dependent and changes across the sky. Therefore it is not necessarily well monitored by blank-field observations taken at different positions on the sky.

A description of the following method to correct for each background component can also be found in Majerowicz & Neumann (2001) and Pratt et al. (2001). A similar method was used by Markevitch & Vikhlinin (2001) for Chandra data.

In the data sets of Abell 1835 and of the blank-sky observations, which we use for background subtraction, we attribute to each event a weighting factor to correct for vignetting (see section 3.2). This is obviously wrong in case of particle events, which are not vignetted. If the particle background (after flare rejection, see section 3.1) is constant with time in terms of spectrum shape, what we assume here, we can correct for the falsely attributed weighting factor by subtracting the background spectrum in the same detector region as the source spectrum. In order to account for possible long time intensity variations in the particle induced background, we determine the count rates of source and background observation at high energies ($`>10`$ keV), where the quantum efficiency for X-ray photons are low. The ratio of the countrates is used to normalize the corresponding background in the spectral fits. For the mos2 we use the 10-12 keV band for determining the background intensity, and for the pn-data the 12-14 keV band, respectively. In our case, a normalization factor of 1 is found. With this method we correct for the particle background. However we do not yet necessarily account for the spatial variations of the cxb.

Abell 1835 is located in a peculiar region of the sky, where an excess of low energy cxb can be seen (Snowden et al. 1997). To correct for this local excess of the cxb, we use regions of the detector in which x–ray emission from the cluster is negligible. We assume in the following that the cxb does not change significantly across the 30′ of detector field of view. We extract a spectrum in the annulus between 7′ and 12′ from the cluster center. From this spectrum outside of the cluster, we subtract the blank-sky spectrum extracted from the same detector region. The residual spectrum represents the difference of the cxb between the observation of Abell 1835 and the blank-sky observations. These residuals were subsequently subtracted from the cluster spectrum, from which we already subtracted the blank-sky spectrum in the same detector regions. We thus yield spectra which are now properly corrected for the two background components described above. In figure 5, we present background corrected spectra from mos2 and pn data which will be used in section 4.4 for the estimate of the mean temperature outside the cooling flow region.

For the pn data there exists another source of contamination, the so called out of time (oot) events, i.e. events counted during the read-out (see Strüder et al. 2001). They are visible as a jet–like structure in figure 4. To estimate the impact of these oot events on the temperature determination, we extract spectra in regions defined in table 1 in which we exclude the oot affected regions. Comparing the fit results of those oot corrected spectra with the non-corrected oot spectra we find a slight increase in the temperature estimates of a few percent, however, always well within the error bars. We thus conclude that the oot’s do not play an important role in the temperature determination for this particular observation. In the following analysis of pn data, we use the non-oot corrected spectra.

#### 4.2.2 Temperature profile

Abell 1835 appears to be a relaxed cluster of galaxies (see section 4.1). Therefore we assume that the temperature structure of this cluster is spherically symmetric. We extract cluster spectra in concentric annuli around the cluster center. The size we choose for the annuli is a compromise between number of source counts to have sufficient statistics on one hand and small enough regions to determine the temperature profile with small bins on the other hand. In our analysis, we exclude the most luminous point sources which are located inside the field of view of the cameras.

We group spectra such that each bin has a signal-to-noise ratio greater than 3$`\sigma `$ after background subtraction. This grouping allows us to assume Gauss statistics, which is important for $`\chi ^2`$ fitting.

For the fit to the spectra, we use xspec version 11.0.1 (Arnaud 1996) and we model the obtained spectra with :

$$\mathrm{Model}=\mathrm{Wabs}(\text{n}\text{H})\times \mathrm{Mekal}(\mathrm{T},\mathrm{A})$$
(1)
where $`\mathrm{Wabs}`$ is a photo electric absorption model (Morrisson & McCammon 1983) and $`\mathrm{Mekal}`$ a single temperature plasma emission model (Mewe et al. 1985 ; Kaastra 1992 ; Liedahl et al. 1995). The free parameters are the galactic hydrogen column density value nH, the x–ray temperature kT and the metallicity (abundance) A.

First we fit all the spectra with a fixed galactic absorption n$`{}_{\mathrm{H}}{}^{}=2.24\times `$1020 cm-2 given by Dickey & Lockman (1990) from 21 cm line measurements. Results are shown in table 1. In a second spectral analysis we leave the nH value as a free fit parameter (see table 2). The resulting fit temperatures agree very well with those obtained with fixed nH which gives us high confidence in our spectral analysis and our background subtraction. Figure 8 shows the fit results for nH which are in good agreement with the measurements from the 21 cm line. Figure 6 shows the resulting temperature profile up to a physical radius of 1.8 Mpc derived from the table 1. Figure 7 shows the metal abundance profile up to 1 Mpc.

### 4.3 Psf corrected temperature profile

Abell 1835 is a medium distant galaxy cluster with a very peaked surface brightness profile (see figure 11) due to the presence of a cooling flow in the centre. Because of this, the point spread function (psf) of Xmm–newton (Griffiths et al. 2002, Ghizzardi 2001) potentially blurs the central regions of the cluster, which might influence the observed temperature profile (see also Markevitch 2002). In order to assess this influence we perform another spectral analysis, in which we take into account the psf of Xmm–newton. We apply the method which is extensively described by Pratt & Arnaud (2002).

For our analysis we assume that the psf is constant across the extent of the cluster and that its shape is energy independent. We calculate the flux redistribution of all temperature bins due to blurring. This computation is based on a convolution of a double $`\beta `$-model, which is fitted to the surface brightness profile and which gave the following parameters (this fit was performed taking the psf into account). Subsequently, based on these model parameters ($`\beta _{\mathrm{in}}=0.49`$, $`r_{c,\mathrm{in}}=0.1`$’, $`\beta _{\mathrm{out}}=0.74`$, $`r_{c,\mathrm{out}}=0.79`$’) we calculate the flux contribution of bin $`i`$ to bin $`j`$. The result can thus be interpreted as a matrix. The corresponding redistribution elements are shown in figure 9 for the pn data. The redistribution for mos2 is similar. As an illustration of the importance of the psf effect in the centre : only 65% of the flux of the second bin ($`i=2`$) are actually observed in the second bin ($`j=2`$). The remaining 35% are distributed in other bins ($`i=2,j2`$).

After the determination of the redistibution or matrix elements we fit simultaneously the different spectra of our six temperature bins. We take into account the different flux contributions of each bin by fixing the ratios of normalization of the fitted spectral models based on our calculated redistribution matrix. Our resulting $`\chi _{\mathrm{red}}^2`$ for the combined fit is $`\chi _{\mathrm{red}}^2=1.05`$ for 4775 degrees of freedom. The resulting temperature profile is displayed in figure 10 as well as in table 3. Briefly, when we take into account the psf we obtain a steeper temperature gradient in the cooling flow region of the cluster. The temperature estimates at larger radii are hardly affected by the psf and are within the error bars of the non-psf corrected analysis. Our study presented here is focussed on the external regions of the cluster, therefore we conclude that the psf does not have an important impact on our analysis (see below).

In the following analysis, we will exclusively use the psf corrected temperature profile.

### 4.4 Mean temperature outside the cooling flow

We use the background subtracted spectra for mos2 and pn data in an annulus between 1′ and 6′ — between 300 and 1800 kpc — around the cluster center to obtain the mean temperature of the cluster outside the cooling flow region. We perform a combined fit using the model described above in equ. (1) with a nH value fixed to 2.24$`\times `$1020 cm-2. We find a temperature of $`7.6\pm 0.4`$ keV and an abundance of 0.25$`\pm 0.05`$ with 90% confidence level.

### 4.5 Surface brightness profile

We extract surface brightness profiles from the cluster data. In order to account for vignetting effects we use again the photon weighting method described in section 3.2. In our specific case using the weight factor for each event to account for vignetting, the count rate $`CR_i`$ in the i-th annulus is defined as :

$$CR_i=\frac{1}{T_{\mathrm{expo}}}\frac{\underset{j=1}{\overset{n}{}}w_j}{A_i}$$
(2)
where $`n`$ is the number of events in the i-th annulus, $`w_j`$ the event weight parameter, $`A_i`$ the surface of the selected annulus and $`T_{\mathrm{expo}}`$ the exposure time. In this case, the Poisson error becomes :

$$\sigma _i=\frac{1}{T_{\mathrm{expo}}}\frac{\sqrt{\underset{j=1}{\overset{n}{}}w_j^2}}{A_i}$$
(3)
Serendipitous sources in the field of view are excluded. Figure 11 shows the surface brightness profile of Abell 1835 derived from the mos2 and pn cameras in the 0.3 to 3 keV energy band.

The surface brightness profile is fitted with a $`\beta `$-model (Cavaliere & Fusco–Femiano 1976) in which the surface brightness $`S`$ is defined as :

$$S(r)=S_0\left(1+\left(\frac{r}{r_c}\right)^2\right)^{3\beta +\frac{1}{2}}+B$$
(4)
where $`S_0`$ is the central intensity, $`r_c`$ the core radius, $`\beta `$ the slope parameter and $`B`$ the background intensity. All these parameters must be fitted.

Abell 1835 hosts a cooling flow in its center (see above). Figure 11 shows the characteristic peak of intensity for cooling flows towards the center (e.g. Fabian 1994). In order to avoid fitting in the cooling flow region we choose to exclude a region with a radius of roughly 0.1 r200 (see equ. (5)) — i.e. 0.7′— around the center for the $`\beta `$-model fit. Indeed Neumann & Arnaud (1999) showed that below 0.1 r200 the surface brightness profiles for nearby clusters show a large dispersion due to cooling flow but they look remarkably similar above 0.1 r200. The best fit parameters are $`\beta =0.704\pm 0.005`$ and a core radius $`r_c=202.3\pm 7.1`$ kpc (i.e. $`r_c=0.681\pm 0.024`$′) for the mos2-data. The uncertainties are 68% confidence level. We also perform a $`\beta `$-model fit with the surface brightness profile extracted from the pn data for consistency. We find in this second case $`\beta =0.708_{0.007}^{+0.012}`$, $`r_c=196.9_{5.3}^{+5.1}`$ kpc (i.e. $`r_c=0.663_{0.018}^{+0.017}`$ ′). These fit parameters are thus in good agreement with those found with the mos2 camera.

Our found parameters of $`\beta `$ and $`r_c`$ are in relatively good agreement with the results found by Neumann & Arnaud (1999) for nearby clusters.

## 5 Mass analysis

We assume in the following that the cluster is virialized up to $`r_{200}`$, the radius in which the mean density of the cluster is 200 times the critical mass density. For the calculation of $`r_{200}=r_{\mathrm{virial}}`$ we use the relation by Evrard et al. (1996) :

$$r_{200}=3690\sqrt{\frac{T}{10\mathrm{keV}}}(1+z)^{\frac{3}{2}}$$
(5)
Using the mean temperature of 7.6 keV (see section 4.4), we find $`r_{200}=2293`$ kpc.

### 5.1 Intra cluster medium mass

Assuming that the radial icm distribution follows a $`\beta `$-model, the electron density $`n_e`$ can be written :

$$n_e(r)=n_{e0}\left(1+\left(\frac{r}{r_c}\right)^2\right)^{\frac{3}{2}\beta }$$
(6)
From the electron density we calculate the total gas density of the icm and determine the icm mass by integrating equ. (6). We obtain a central electron density $`n_{e0}=1.47\pm 0.02\times 10^2`$ cm-3 (68% confidence level) assuming a temperature of 7.6 keV and using the $`\beta `$-model parameters determined in section 4.5. Temperature variations in the uncertainties do not change our values for the electron density. The error bars of $`n_{e0}`$ are calculated from the errors of the $`\beta `$-model parameters. We verify that the cluster temperature has no effect on the estimate of $`n_{e0}`$ with the selected energy range.

The icm mass profile is shown in figure 12. We find M$`{}_{\mathrm{gas}}{}^{}=2.56\pm 0.08\times 10^{14}`$ M within $`r=1.7`$ Mpc or 0.75 $`r_{200}`$ with a confidence level of 68%.

We want to stress that we underestimate the central gas density in the centre due to neglecting the cooling flow excess. However, this underestimation at small radii does not play a significant role when we integrate the gas mass out to large radii.

In fact, when we fit the surface brightness profile (figure 11) so that we only account for the inner 210 kpc, we find $`\beta `$=0.49, $`r_c`$=52 kpc and $`n_{e0}=5.55\times 10^2`$ cm-3. By integrating equ. (6) with these parameters, we find a gas mass at r=0.1 $`r_{200}`$, which is 18% higher than the one using the global $`\beta `$-model fit parameters. At $`r=0.75r_{200}`$, which is the radius up to which we can constrain the temperature and thus total mass profile we underestimate the gas mass by 1% due to neglecting the gas mass linked to the central gas excess with respect to our global $`\beta `$-model , which is negligible.

### 5.2 Cluster total mass

Assuming spherical symmetry and hydrostatic equilibrium, we calculate the gravitational mass of the cluster Abell 1835 with the hydrostatic equation :

$$\mathrm{M}_{\mathrm{tot}}(<r)=\frac{k\mathrm{T}}{G\mu m_p}r\left(\frac{d\mathrm{ln}n_e}{d\mathrm{ln}r}+\frac{d\mathrm{ln}\mathrm{T}}{d\mathrm{ln}r}\right)$$
(7)
where $`k`$ is the Boltzmann constant, T the gas temperature, $`G`$ the gravitational constant, $`\mu `$ the mean molecular weight of the gas ($`\mu 0.6`$), $`m_p`$ the proton mass and $`n_e`$ the electronic density. Including equ. (6) in equ. (7), M can be written :

$$\mathrm{M}_{\mathrm{tot}}(<r)=\frac{k}{G\mu m_p}r^2\left(\frac{d\mathrm{T}}{dr}3\beta \mathrm{T}\frac{r}{r^2+r_c^2}\right)$$
(8)
In order to account for the error bars in the temperature profile, we use a Monte–Carlo method described by Neumann & Böhringer (1995). This method allows to transform error bars of the temperature profile into error bars of the mass profile. This transformation is done by calculating randomly temperature profiles, which fit in the actual measured temperature profile and to subsequently calculate the mass profile corresponding to the randomly determined temperature profile. We determine 10000 random temperature profiles. From the corresponding 10000 mass profiles we calculate the mean mass profile and the errors. The uncertainties are determined by taking into account the uncertainties of the temperature profile (90% confidence level) and the standard deviation of the mass at a given radius. The resulting errors are calculated so that they correspond to 1$`\sigma `$ error bars (68% confidence level). The mass profile is shown in Fig.13. At a radius corresponding to 75% of the virial radius (or 1.7 Mpc) we obtain a total mass333For a comparison, the total mass is $`1.14\pm 0.14\times 10^{15}`$ M when we use the temperature profile without psf correction (table 1). of : $`M_{tot}=1.23\pm 0.18\times 10^{15}`$ M (68% confidence level).

Assuming that Abell 1835 is isothermal with a temperature of 7.6 keV (see section 4.4), equ. (8) can be written as :

$$\mathrm{M}_{\mathrm{tot}}(<r)=\frac{3k\beta }{G\mu m_p}\mathrm{T}\frac{r^3}{r^2+r_c^2}$$
(9)
In this case taking into account error bars on T, $`\beta `$ and $`r_c`$ parameters, the total mass of Abell 1835 is $`1.00\pm 0.06\times 10^{15}`$ M which is consistent with the mass obtained using the Monte-Carlo approach. This was expected since the temperature profile of Abell 1835 (see figure 10) seems to be isothermal at high radii.

### 5.3 Gas mass fraction

The gas mass fraction is simply the ratio of the icm mass to the total mass. The profile of the gas mass fraction of Abell 1835 is shown in figure 14. We do not display the gas mass fraction below 10% of the virial radius because of the cooling flow region and our underestimate of the gas mass there (see section 4.5). The errors bar are calculated with classical error propagation and are at a confidence level of 68%. Within $`r=0.75r_{200}`$, we find $`f_{gas}h_{50}^{1.5}=0.207\pm 0.037`$ where $`h_{50}`$ is the present time Hubble constant in units of 50 km/s/Mpc.

## 6 Discussion

### 6.1 The gas mass fraction and implications for cosmology

Abell 1835 is an ideal cluster to determine the gas mass profile with the hydrostatic approach, since the Xmm–newton image does not show indications of substructure at large scales444Schmidt et al. (2001) reported small scale substructure in the Chandra data in the cooling flow region, which does not influence our results on large scales., which might contaminate the mass determination.

Our determined gas mass fraction profile of Abell 1835, which we determine up to 0.75 $`r_{200}`$ is essentially flat. The mean value is about 21% and within 0.75 $`r_{200}`$ we obtain a gas mass fraction of $`0.207\pm 0.037h_{50}^{3/2}`$. We can take this gas mass fraction as lower limit for the baryon fraction in the cluster because we do not take into account other contributions of baryons, such as stars or interstellar medium. Extrapolating this lower limit through out the universe (which was shown to be valid by White et al. (1993) and references therein), we can determine the cosmological matter density parameter $`\mathrm{\Omega }_m`$ by combining our result with studies based on primordial nucleosynthesis, which give firm constraints on the density parameter of baryons in the universe. Recent results give $`\mathrm{\Omega }_bh_{50}^2=0.080\pm 0.008`$ (Burles et al. 2001). We can thus apply :

$$\frac{M_{\mathrm{gas}}}{M_{\mathrm{tot}}}=0.207\pm 0.037h_{50}^{3/2}=\frac{\mathrm{\Omega }_\mathrm{b}}{\mathrm{\Omega }_\mathrm{m}}=\frac{0.08\pm 0.008h_{50}^2}{\mathrm{\Omega }_\mathrm{m}}$$
(10)
If the icm was the only source of baryons in clusters, above equation gives us $`\mathrm{\Omega }_\mathrm{m}=0.386\pm 0.108`$ (68% confidence level). However, since there exist other baryon sources in Abell 1835 such as galaxies, this implies that our result is an upper limit with : $`\mathrm{\Omega }_\mathrm{m}<0.5h_{50}^{1/2}`$. This result is in very good agreement with other recent results based on cmb measurements combined with redshift surveys (Efstathiou et al. 2001) and studies on distant supernovæ (Perlmutter et al. 1999) which favor a low density universe with a cosmological constant. If we additionally knew the baryon fraction in galaxies up to high precision, we could put even tighter constraints on $`\mathrm{\Omega }_m`$.

### 6.2 The temperature profile – comparison with other recent studies

A lot of effort has been spent already for the determination of temperature profiles in the icm. The corresponding results vary from each study to another quite significantly. There have been reported declining temperature profiles, with a more or less universal shape based on Asca and Beppo-sax data (Markevitch et al. 1998, de Grandi & Molendi 2001). Other studies based on data from the same instruments, on the contrary, find flat temperature profiles or profiles with large dispersions (Irwin & Bregman 2000, White et al 2000).

All these studies are limited to radii smaller than 0.5 $`r_{200}`$ due to instrumental limitations. Furthermore, studies based on Asca and Beppo-sax were subject to possible systematic uncertainties due to the large point spread function of the mirrors, which additionally can vary in size depending on energy. Chandra results presented by Allen, Schmidt & Fabian (2001) suggest that massive, relaxed clusters are isothermal up to $`r_{2500}`$ and show a universal temperature decline in the center linked to the cooling flow.

Xmm–newton with its better spatial resolution than Asca and Beppo-sax and larger effective area than comparable observatories allows to a) better constrain the temperature profiles and b) to trace them out to larger radii. So far studies based on Xmm–newton show either flat profiles (Arnaud et al. 2001a ; Tamura et al. 2001) for hot clusters, which is in agreement with the results by Allen et al. (2001) on Chandra data of clusters. A Xmm–newton study on the cool cluster Sérsic 159-03, however, shows a dramatic temperature decline in the outer parts of the cluster (Kaastra et al. 2001).

Our results suggest strongly a flat temperature profile with a decline towards the center linked to the cooling flow present in Abell 1835 (Peterson et al. 2001, Schmidt et al. 2001 and references therein). This profile fits well with the empirically found model by Allen et al. (2001), which is valid up to $`r_{2500}`$, which corresponds roughly to $`0.3\times r_{200}`$ and which is displayed in Fig 15. However, an individual study undertaken by Schmidt et al. (2001) on Abell 1835 observed with Chandra shows a temperature profile which is increasing with radius. At radii larger than 400 kpc Schmidt et al. (2001) find temperature estimates of about 12 keV, which is inconsistent with our Xmm–newton results. The fact that we find good agreement between the individual spectral modeling of mos2 and pn cameras gives us confidence in our results.

While undertaking our study concerning the background subtraction we observed that minor differences in the extraction of flares between background and source observation can lead to declining or increasing temperature profiles. Doing monitoring on the background level we find that quiescent periods do not always have the same background level and have sometimes, even though rarely, intensities close or above the threshold criterion to extract flares. In order to screen the data properly for background variations and to determine the level of the particle (instrumental) background accurately it is important to look at parts of the observation where no astrophysical source emission can be detected like the high energy part of the spectrum, where only particle background can be seen (the particle background has a hard spectrum and the effective area of Xmm–newton at energies above 10 keV for x–ray photons is relatively small). It is important to compare and to apply exactly the same flare rejection for both the source and the background observation to avoid wrong background subtraction, which can bias the temperature estimates. Observing a temperature profile which increases with radius can be an indication of a slightly higher particle background in the source observation than in the background observation.

The presence of remaining solar flares in the analysis of Schmidt et al. (2001) was recently discussed by Markevitch (2002), who reanalyzed the Chandra data of Abell 1835. Markevitch (2002) found that a temperature profile analysis based on a more thorough flare treatment leads to a basically flat temperature profile in the outer parts of the cluster, which agrees well with our derived temperature profile (see figure 16) at large scales.

### 6.3 Global temperature estimates

There is some dispersion observable between past temperature determinations of Abell 1835 based on Asca-data. Allen & Fabian (1998) found a temperature of $`9.8_{1.3}^{+2.3}`$ keV, based on Asca-sis and gis-data, including a multiphase cooling flow model, which is now known not to be an accurate description for the cooling flow of Abell 1835 (Peterson et al. 2001). White (2000) only used Asca-gis data and obtained an overall temperature of kT=$`8.5_{0.5}^{+1.5}`$ keV when applying a multiphase cooling flow model and kT=$`8.2\pm 0.5`$ keV when using a single temperature model. Both results are consistent within the error bars with our results of kT=7.6$`\pm 0.4`$ keV. Furthermore, Mushotzky & Scharf (1997) found a temperature estimate of kT=$`8.2\pm 0.5`$ keV, which is again comparable with our results and the ones of White (2000). The results of Allen & Fabian (1998) agree within the error bars with other results based on Asca, but are higher than our mean temperature. While our temperature results are somewhat lower than the results of Asca, the Chandra results for the temperature of Abell 1835 by Schmidt et al. (2001) are above the results based on Asca.

## 7 Conclusion

Xmm–newton data allow us to measure the temperature profile of Abell 1835 up to 0.75 $`r_{200}`$. In order to determine the temperature profile of the cluster accurately we apply a method, which corrects for the various kinds of background and a photon weighting method which allows to correct for vignetting effects. We correct the temperature profile for psf-effects and see that the resulting profile is not affected by the psf at large radii.

Our resulting psf corrected temperature profile is consistent with being flat in the outer parts. The mean cluster temperature is 7.6 keV outside the cooling flow region. In the central parts, below a radius of 1′(300 kpc) we see a temperature decline linked to the cooling flow of the cluster. We fit a $`\beta `$-model to the outer radii of the surface brightness profile of Abell 1835 and obtain the best fit results of $`\beta =0.704\pm 0.005`$ and $`r_c=202.3\pm 7.1`$ kpc. We use these best fit parameters as input for the hydrostatic approach to calculate the mass profile of the cluster and find a mass at 0.75 $`r_{200}`$ of $`1.23\pm 0.18\times 10^{15}`$ M. The corresponding gas mass fraction, which is constant with radius at radii $`r>0.1r_{200}`$, is with 20.7%$`\pm 3.7`$ comparable with the results found in other studies.

Our global temperature estimate for Abell 1835 based on Xmm–newton is somewhat low when compared to other results, but agree in the error bars with other studies based on Asca-data. We assume that these discrepencies are linked to small remaining uncertainties in the calibration of Xmm–newton , and on which is still being worked on.

###### Acknowledgements.

We would like to thank Monique Arnaud for useful discussions, suggestions on the background subtraction, and a very important participation in the development of the codes we use here. Equally we are greatful to René Gastaud, who wrote most of the codes concerning the vignetting correction. Furthermore we would like to thank D. Lumb for his effort to accumulate blank-sky observations, which is essential for our study presented here. We also thank the referee for useful remarks. T. H. R. thanks the Celerity Foundation for support. T. H. R. was also supported by NASA XMM-Newton Grant NAG5-10075. Also, we are very grateful to the entire epic-calibration team, which shared all essential information with us, which made this study possible. The results presented here are based on observations obtained with Xmm–newton, an Esa science mission with intruments and contributions directly funded by Esa Member States and the Usa (Nasa).


format: text

text:
# The dynamical foundation of fractal stream chemistry: The origin of extremely long retention times

## I Introduction

The distribution of travel times required for chemical substances to travel through a catchment is an important hydraulic quantity, determining the retention of pollutants until they are eventually released to a stream or lake. The same distribution controls the transport of chemicals in subsurface hydrological systems. The nature of this distribution determines the expected ecological impact of contaminants.

In a long-term catchment study conducted at Plynlimon, Wales, Kirchner et al.  (hereafter denoted KFN) compare the (input) time series of the concentration $`c_R(t)`$ of chloride tracer in the rainfall to the (output) time series of the concentration $`c_s(t)`$ into the Hafren stream. They relate the concentrations through the convolution integral

$$c_s(t)=_0^{\mathrm{}}h(t^{})c_R(tt^{})𝑑t^{}$$
(1)
where the effective travel time distribution $`h(t)`$ governs the lag time between injection of the tracer through rainfall and outflow to the stream.

KFN observe that the streamflow (volume of water) responds rapidly to storm rainfall inputs while the $`c_s(t)`$ response is highly damped (i.e., low-level contaminant delivery to streams for a very long time). The timescales of catchment hydrologic and chemical response are resolved using spectral methods Duffy and Gelhar, (1985) (in which the input/output response can be compared at each (time) wavelength). The spectral power of the water fluxes of rainfall and streamflow coincide while the spectra of $`c_s(t)`$ show strong attenuation on all wavelengths less than 5-10 years. Further, the $`c_R(t)`$ spectra scale as white noise in sharp contrast to the fractal $`1/f`$ scaling of the $`c_s(t)`$ spectra. Using these results in the Fourier transform of the convolution in (1), with $`f`$ denoting frequency, and $`C_s`$, $`H`$ and $`C_R`$ representing the respective Fourier transforms,

$$C_s(f)=H(f)C_R(f)$$
(2)
and the relation of the power spectra is

$$|C_s(f)|^2=|H(f)|^2|C_R(f)|^2$$
(3)
Since $`|C_R(f)|^2`$ is nearly a white noise spectrum (a constant) one has $`|H(f)|^2|C_s(f)|^2`$. From the measured power-law scaling of $`|C_s(f)|^2f^{0.97}`$, KFN conclude that

$$h(t)t^m$$
(4)
(from $`H(f)f^{(1m)}`$), where $`m0.5`$.

In order to ensure that $`h(t)`$ is integrable and possesses a finite average travel time, KFN arbitrarily chose the gamma distribution

$$h(t)=\frac{t^{\gamma 1}}{\varrho ^\gamma \mathrm{\Gamma }(\gamma )}e^{t/\varrho }$$
(5)
with $`\gamma =0.48`$ and $`\varrho =1.9`$ years, with a characteristic travel time $`T_h_0^{\mathrm{}}h(t)t𝑑t=\gamma \varrho `$. Due to the relatively large value of $`\varrho `$, essentially the entire data correspond to the power-law $`h(t)t^{\gamma 1}`$ with the power spectrum $`|C_s(f)|^2|H(f)|^2f^{2\gamma }`$. KFN refer to similar scaling between $`f^{0.7}`$ and $`f^{1.2}`$ found in Scandinavian and North American field sites, which indicates a certain ubiquity to fractal, scale-free forms (which is the issue we address in this letter).

## II Transport models

The first step to understanding transport in the catchment is the clarification of the meaning of (1). The $`h(t)`$ is the effective response to a pulse of rain falling on the entire area of the catchment. Every point of this area is a source of chloride and the stream is a line sink for the chloride. We simplify the area to be a rectangle of width $`\lambda `$ about this stream sink (with periodic boundary conditions on the sides perpendicular to the streamScher and Montroll, (1975)). The sink or absorbing boundary ensures the correct “counting rate” of chloride at the arriving point, i.e., the first passage time distribution $`F(t;l)`$ is the appropriate travel time distribution from a pulse source at point $`l`$. In terms of this intrinsic distribution,

$$c_s(t;l_s)=_0^{\mathrm{}}\underset{l\mathrm{\Omega }}{}F(t^{};ll_s)c_R(tt^{};l)dt^{}$$
(6)
where $`c_s(t;l_s)`$ is the chloride concentration at the stream position $`l_s`$ at time $`t`$, $`\mathrm{\Omega }`$ is the size of the catchment and $`c_R(t;l)`$ the rain-input at a position $`l`$ in $`\mathrm{\Omega }`$. We can consider the sampling positions $`\{l_s\}`$ of $`c_s`$ to be a small region compared to $`\mathrm{\Omega }`$ and hence $`c_s(t;l_s)c_s(t)`$. To do the $`l`$-sum we assume that the rainfall is distributed uniformly in $`\mathrm{\Omega }`$, $`c_R(t;l)c_R(t)`$. Hence, we recover (1) now with

$$h(t)\underset{l\mathrm{\Omega }}{}F(t;l)$$
(7)
The distribution $`F(t;l)`$ must be integrable in time and thus so $`h(t)`$ for a finite $`\mathrm{\Omega }`$. The basis of comparison for various transport approaches is the computation in (7).

### II.1 Advection-dispersion equation (ADE)

The textbook approach to transport of passive tracers in both surface and subsurface systems usually focuses on the advection-dispersion equation (ADE)

$$𝒲/t+v𝒲/x=D^2𝒲/x^2$$
(8)
e.g., for 1d with constant $`v`$ the average fluid velocity and $`D`$ the dispersion coefficient. This equation governs the temporal evolution of the probability density function (pdf) $`𝒲(x,t)`$ of finding a tracer particle, which undergoes dispersive motion, at a certain position $`x`$ at time $`t`$ after release at $`t=0`$. The pdf $`𝒲`$ is the normalized concentration profile.

The travel time distribution belonging to (8) can be obtained by the Laplace Transform ($`LT`$) technique. This yields the $`LT`$ $`\stackrel{~}{}(u;x)_0^{\mathrm{}}(t;x)e^{ut}𝑑t`$ of the first passage time from $`x`$ to the absorbing boundary at $`x=0`$, $`\stackrel{~}{}(u;x)=\mathrm{exp}\left(x\left[v\sqrt{4Du+v^2}\right]/(2D)\right)`$. $`LT`$ inversion delivers the travel time distribution

$$(t;x)=\frac{x}{\sqrt{4\pi D}t^{\frac{3}{2}}}\mathrm{exp}\left(\frac{(vtx)^2}{4Dt}\right)$$
(9)
which when inserted for $`F(t;x)`$ in (7) yields

$`h(t)`$ $``$ $`{\displaystyle _0^\lambda }(t;x)dx={\displaystyle \frac{v}{2}}[\mathrm{erf}(z)+\mathrm{erf}({\displaystyle \frac{P_e}{4z}}z)`$ (10)
$`+`$ $`{\displaystyle \frac{\mathrm{exp}(z^2)\mathrm{exp}\left((\frac{P_e}{4z}z)^2\right)}{\sqrt{\pi }z}}]`$
where the Peclet number $`P_e\lambda v/D`$, $`z\frac{1}{2}\sqrt{v^2t/D}`$ and erf$`(z)`$ is the error function (the same result (10) was also derived in Kirchner et al., 2001). For $`t<0.1D/v^2`$, $`h(t)t^{\frac{1}{2}}`$ which is similar to (4) for $`m0.5`$. \[This upper limit on the time range holds for $`P_e1`$; for $`P_e<1`$ the limit is smaller.\] The range is estimated by assuming $`D`$ to be the macrodispersion, $`D=v\alpha `$, where $`\alpha `$ is the dispersivity. Assuming $`v10^2`$ m/yr and on the km scale $`\alpha (110^2)`$ m, $`t<0.1\alpha /v(10^310^1)`$ yr. One needs $`v10`$ m/yr and $`\alpha 10^2`$ m in (10), for $`h(t)t^{\frac{1}{2}}`$ to partially overlap the time range of observation of this dependence Kirchner et al., (2000). Moreover, this behavior is only the same as (4) in the special case where $`m`$ (or $`\gamma `$) = $`\frac{1}{2}`$. The catchment data (spectral power $`f^{2\gamma }`$) quoted above cover a range of $`0.72\gamma 1.2`$.

For large $`t`$, $`h(t)`$ in (10) exhibits an exponential decrease $`t^{\frac{3}{2}}`$ $`\mathrm{exp}\left(v^2(4D)^1t\right)`$ which is faster than (5). The decrease in $`h(t)`$ at large $`t`$ ($`t>>\lambda /v`$) is due to the finite size of the catchment and not an arbitrary limiting time. The relation between the operative time range of (4) and $`\lambda `$ will be determined in the following.

### II.2 Continuous-time random walk and time-fractional-ADE

As shown above the standard treatment of the transport using the ADE yields a marginal accounting of the data. A transport model must determine a $`F(t;l)`$ that when inserted in (7) gives rise to (4) for over three decades of time Kirchner et al., (2000). As we will show this is best expedited with a non-Gaussian form of $`F(t;l)`$. A particularly dispersive form of $`F(t;l)`$ is associated with anomalous transport which has been observed in a wide range of disordered and complex systems: electron hopping/multiple-trapping in amorphous semiconductors Scher and Montroll, (1975), particle migration in fracture networks Kosakowski et al., (2001), contaminant transport in heterogeneous porous media Berkowitz and Scher, (1998), anomalous diffusionMetzler and Klafter, (2000), and Hamiltonian chaosKlafter et al., (1996). There are key common features in these anomalous transport phenomena, e.g., non-Gaussian propagation of an initial pulse of particles with a mean position $`\mathrm{}(t)`$ and standard deviation $`\sigma (t)`$ exhibiting the same sublinear dependence on $`t`$ (in the presence of a bias or head). A continuous time random walk (CTRW) transport process has successfully accounted for these unusual basic features. The anomalous behavior is an outcome of a CTRW governed by a long-tail distribution of the individual transition times or event times (which have to be defined in each physical context, e.g., trap release),

$$\psi (t)A_\beta t^{1\beta }$$
(11)
where $`A_\beta `$ is a constant. Over the observation time range in which (11) obtains we have for the propagating plume $`P(l,t)`$

$$\mathrm{}(t)t^\beta ,\sigma (t)t^\beta 0<\beta <1$$
(12)
The highly dispersive nature of this $`P(l,t)`$ can be discerned by the ratio $`\sigma (t)/\mathrm{}(t)`$ constant, instead of the familiar $`1/\sqrt{t}`$ for a Gaussian plume. The latter is produced for $`\beta >2`$ (while for $`1<\beta <2`$ one has an intermediate $`t`$-dependence).

The CTRW accounts naturally for the cumulative effects of a sequence of transitions which constitutes the particle transport. The formalism of the CTRW has been well documented in the literature \[Berkowitz and Scher, 1998; Metzler and Klafter, 2000\]. For brevity we show the key equation for our purpose here in Laplace space (for three dimensions)

$$\stackrel{~}{F}(𝐬,u)=\stackrel{~}{P}(𝐬,u)/\stackrel{~}{P}(0,u)𝐬0$$
(13)
where the particle starts at the origin at $`t=0^+`$, $`\stackrel{~}{P}(𝐬,u)`$ is the $`LT`$ of $`P(𝐬,t)`$, the probability density to find the particle on $`𝐬`$ at time $`t`$ (plume), and $`\stackrel{~}{F}(𝐬,u)`$ is the $`LT`$ of $`F(𝐬,t)`$, the probability per time for the particle to first arrive at $`𝐬`$ at time $`t`$.

Explicit evaluations for $`\stackrel{~}{F}(𝐬,u)`$ have been determined for symmetric situations (e.g., plane (line) source to plane (line) sink) and with the use of (11) Scher and Montroll, (1975). The symmetric cases allow a reduction of the problem to a one-dimensional form of $`F`$ which we write as $`F(t;x)`$ and

$$\stackrel{~}{F}(u;x)=\mathrm{exp}(xu^\beta /\overline{l})(0<\beta <1)$$
(14)
($`u`$ dimensionless). The inverse $`LT`$ of $`\stackrel{~}{F}(u;x)`$ has been expressed in terms of a class of Fox $`H`$-functions Metzler and Klafter, (2000), but it is more expedient to work directly with (14) in the evaluation of $`h(t)`$, which will be considered in the next section.

Mathematically, the CTRW (in the special case using (11) and $`(0<\beta <1)`$) is interchangeable with the time-fractional advection-dispersion equation (FADE) Metzler and Klafter, (2000)

$$\frac{W}{t}=_0D_t^{1\beta }\left(v_\beta \frac{}{x}+d_\beta \frac{^2}{x^2}\right)W(x,t)$$
(15)
where the Riemann-Liouville operator is defined in terms of the convolution Metzler and Klafter, (2000),

$$_0D_t^{1\beta }W(x,t)=\frac{1}{\mathrm{\Gamma }(\beta )}\frac{d}{dt}_0^tdt^{}\frac{W(x,t^{})}{(tt^{})^{1\beta }},\mathrm{\hspace{0.17em}0}<\beta <1.$$
(16)
The form of FADE (15) is the natural generalization of the ADE (8) for power-law forms of $`\psi `$. The solutions of FADE reproduce those of the CTRW in this special case, and in the spatial continuum limit they form a natural bridge between the CTRW and the ADE Berkowitz and Scher, (2001).

## III Long Retention Times

We proceed to insert (14) into the $`LT`$ of (7) to obtain

$$\stackrel{~}{h}(u)_0^\lambda \mathrm{exp}(lu^\beta /\overline{l})𝑑l=\overline{l}u^\beta (1\mathrm{exp}(\lambda u^\beta /\overline{l}))$$
(17)
where $`\overline{l}`$ is the mean step distance. The expression for $`\stackrel{~}{h}(u)`$ in (17) has been thoroughly studied in another context Scher and Montroll, (1975); the main features are

$$h(\tau )\{\begin{array}{cc}\tau ^{\beta 1},\hfill & \tau <\tau ^{}\hfill \\ \tau ^{\beta 1},\hfill & \tau >\tau ^{}\hfill \end{array}$$
(18)
($`\tau `$ dimensionless time). The exponent for $`\tau >\tau ^{}`$ ensures that $`h(\tau )`$ is integrable. On a log-log plot $`h(\tau )`$ is two constant slopes, $`\beta 1`$, $`1\beta `$, with a turnover range between them. The center time of this range, $`\tau ^{}`$, can be estimated as the time for the argument of the exponent in (17) to be $`O(1)`$ (using $`u1/\tau `$)

$$\tau ^{}(\frac{1\beta }{\beta })^{\frac{1\beta }{\beta }}(\beta \lambda /\overline{l})^{\frac{1}{\beta }}$$
(19)
(The $`\beta `$-factors derive from a more detailed analysis \[Scher and Montroll, 1975, Appendix C\].) For $`\beta =\frac{1}{2}`$, in (17) one can do the inverse $`LT`$ and obtain $`h(\tau )(\pi \tau )^{\frac{1}{2}}(1\mathrm{exp}(\lambda ^2/4\overline{l}^2\tau ))`$, clearly showing (18)-(19).

The change from $`\tau `$ to $`t`$ is model dependent (to be discussed below). We use $`\tau =v_ot/\overline{l}\epsilon `$, where $`v_o`$ is a characteristic velocity (of the velocity distribution) and $`\epsilon `$ is a constant. For $`\beta =\frac{1}{2}`$, $`v_o`$ 100 m/yr, $`\lambda \frac{1}{2}`$ km, $`\overline{l}`$ 30 m and $`\epsilon \frac{1}{2}`$, one has (from (19) with $`\tau ^{}=v_ot^{}/\overline{l}\epsilon `$), $`t^{}`$ 10 yr, which is a reasonable time scale for the change from $`h(t)t^{\beta 1}`$. Hence, our scaling result (18) agrees with the KFN data (4) over the measurement time range ($`>`$ 3 decades for t$``$ 10 yrs) with $`\beta =\gamma `$ ($`1m`$). The $`t`$-“cutoff” for $`h(t)`$ in (18) is algebraic, $`t^{\beta 1}`$, not exponential as in (5)! This behavior of $`h(t)`$ is indicative of extremely long chemical retention times in catchments. The turnover to the $`t^{\beta 1}`$ dependence (in (18)) is a prediction of the CTRW theory and has not yet been observed. A finite value of $`T_h`$ depends on a change in the behavior of (11). For $`tt^{}`$, the tail of $`\psi (t)`$ becomes steeper and for $`\beta >2`$, the intrinsic arrival time distribution evolves to $`(t;x)`$.

Independent tests of our transport model can be obtained with experiments using a different tracer (d) and boundary conditions. For injection of the tracer d at various catchment points $`\{𝐱\}`$, the concentration in the stream $`c_d(t;𝐥_𝐬)`$ is proportional to $`F(𝐱𝐥_𝐬,t)`$, determined for these boundary conditions. \[This result answers the question raised by Stark and Stieglitz  about a site-specific spill.\] This is also a method for obtaining a value of $`\beta `$, which determines the shape of $`F(𝐱,t)`$.

There are many different physical mechanisms (cf. Refs. above) that can give rise to the $`\psi (t)`$ in (11) that has generated our results for $`h(t)`$. The common features are the representation of the transport process as a series of transitions \[Berkowitz and Scher, 1998, 2001\] and the encounter within this series of a sufficient number of transition-times that are much larger than the median one. These relatively few long-time events, which can be due to release from traps/stagnant-regions and/or passage through low velocity segments, can have a dominant influence on the overall transport. The relative weighting of these events is expressed by the exponent of the power-law in (11).

The hydrologic environments beneath the catchment subregions can contain a sufficient number of these low velocity sections and stagnant dead-ends. The local slope (and interconnection Rodriguez-Iturbe and Rinaldo, (1997)) of the catchment side-channels also affects the velocity distribution. The subsurface of the catchment basin can be modeled as a heterogeneous porous medium and/or a random fracture network Berkowitz and Scher, (1998). In the latter the fragment-length distribution, $`f(s)`$, and the fragment-velocity histogram, $`\mathrm{\Phi }(\xi )`$, are mapped onto the probability rate for a transit time t (through a fragment) with a displacement s,

$$\psi (𝐬,t)\mathrm{\Phi }(\xi )f(s)$$
(20)
where $`|\xi |=1/v`$, $`\widehat{\xi }=\widehat{𝐯}`$, $`t=s\xi `$ and $`\psi (t)_s\psi (𝐬,t)`$. The anomalous transport observed in the fracture network particle tracking simulations is due to the power-law tail of $`\mathrm{\Phi }(\xi )\xi ^{1\beta }`$ at large $`\xi `$. This is the expected behavior to be found in the catchment subsurface $`\mathrm{\Phi }_c(\xi )`$ distribution because this hydrologic environment is the same as ones modeled successfully by the CTRW theory \[Berkowitz and Scher, 1998, 2001\]. The experiments outlined above can test these properties.

## IV Conclusion

In this letter we introduce the phenomenon of anomalous transport into an analysis of the comparison of the power spectra of the time series of chloride in rainfall and in (stream) runoff. An effective distribution of travel times $`h(t)`$ is derived (by KFN from the data) which we show is composed of the summation from all sites of the source distribution $`F(𝐬,t)`$ in the catchment. The latter, in turn, is composed of a sequence of transitions governed by a power-law (tail) distribution (11). The results agree naturally with the data reported from a number of catchment studies, in the sense of a power-law scaling of $`h(t)`$ over decades of $`t`$ in the observational time range (using reasonable parameters). A consequence of our transport model is a power-law ”cutoff” of $`h(t)`$ predicting extremely long chemical retention times with subsequent impact on the long-term effects of contaminants on these ecosystems.

###### Acknowledgements.

We thank James W. Kirchner and Colin Stark for a critical reading of the manuscript. RM acknowledges the DFG and the European Commission within the Emmy Noether and SISITOMAS programmes, and BB acknowledges the European Commission (Contract EVK1-CT-2000-00062), for financial support.


text:
# Series expansion of the quark determinant in the number of quark-antiquark pairs

## 1 Introduction

A clearcut separation of the valence and sea quark contributions to the QCD partition function might help in the study of a number of issues, starting with the spin content of the nucleon  up to QCD at finite baryonic density . The identification of the different contributions can be achieved by means of an expansion of the partition function with respect to the number of quark-antiquark pairs, which is the goal of the present paper. But the practical usefulness of such an expansion depends on the number of terms necessary to get the desired accuracy.

In the case of the nucleon, for instance, only a few quarks of the sea are expected to contribute significantly. If the partition function is split in the corresponding terms only a few will be important.

At first sight the case of QCD at nonzero baryon density appears quite different, because in principle one is interested in the termodynamic limit. To investigate this limit an infinite number of terms would obviously be necessary, making our expansion of little use. In actual numerical simulations, however, the physical volumes accessible are of the order of a few $`fm^3`$. If we consider that normal nuclear density is about $`3quarks/6fm^3`$, with 6 valence quarks we are already at about twice the ordinary nuclear density, and it is reasonable to expect that the number of important quark-antiquark pairs be comparable to the number of valence quarks at least at moderate density and temperature.

To our knowledge there exists no clearcut separation of the contributions we are discussing. To tackle this problem we need, to start with, a formulation of the partition function in a given baryonic sector. One such formulation has been constructed by introducing an imaginary chemical potential in the grand canonical partition function and by projecting in a given baryonic sector by a Fourier transformation . But the result is not very transparent with respect to the distinction of the different contributions we are interested in.

We consider here a different form of projection, which leads to a series expansion for the quark determinant. The terms of this expansion have the desired physical meaning: The r-th term gives the gauge-invariant contribution of the valence quarks plus r quark-antiquark pairs. In view of the above considerations, we hope that our expansion will have a reasonable convergence in present numerical simulations.

The technique we employ is fairly general, requiring only the knowledge of the transfer matrix in the Fock space of the quarks. But in this paper we restrict ourselves to the case of Wilson fermions. We hope to be able to treat the Kogut-Susskind fermions in a future work.

We outline the derivation of the series expansion in Sec.2 focussing on the strategy, leaving out many details of the calculations for Sec.3. In Sec.4 we give our conclusions.

## 2 The series expansion of the quark determinant in a given baryonic sector

In this Section we outline the steps leading to a series expansion of the quark determinant in a given baryonic sector. As already said the presentation will be somewhat schematic.

The quark determinant in the absence of any condition on the baryon number is

$$detQ=[d\overline{\psi }\psi ]\mathrm{exp}S_F.$$
(1)
$`S_F`$ is the quark action and $`Q`$ the standard quark matrix which will be spelled out later. Our strategy is to write $`detQ`$ as the trace of the transfer matrix acting in the quark Fock space, impose the restriction to a given baryonic sector, and then rewrite the trace as the determinant of a modified quark matrix. The round trip is done by mapping the Grassmann algebra generated by the quark fields into the Fock space following the construction of Lus̈cher . But while his paper is based on the mere existence of the map, for us it is essential a concrete realization by means of coherent states.

In the transfer matrix formalism the (euclidean) time is treated differently from the spatial coordinates. So we must also treat it differently, and we assume a convention of summation over intrinsic quantum numbers and spatial sites only. These have coordinates $`𝐧`$, while the time will be denoted by $`n_0`$. The link variables $`U_\mu `$ are matrices in color and have spatial matrix elements

$`\left(U_0(n_0)\right)_{𝐦,𝐧}`$ $`=`$ $`U_0(𝐦,m_0)\delta _{𝐦,𝐧}`$
$`\left(U_j(m_0)\right)_{𝐦,𝐧}`$ $`=`$ $`U_j(𝐦,m_0)\delta _{𝐦+𝐣,𝐧}.`$ (2)
The quark field $`\psi `$ carries Dirac, color and flavor indices, denoted altogether by $`\alpha `$, and space-time label $`n`$. We will use with the same convention Grassmann variables $`x`$’s and $`y`$’s and quark creation and annihilation operators $`\widehat{x}^+,\widehat{x},\widehat{y}^+,\widehat{y}`$. But unlike Grassmann variables, the creation and annihilation operators do not carry a temporal index. Thus according to our convention

$`\overline{\psi }(n_0)\psi (n_0)`$ $`=`$ $`{\displaystyle \underset{𝐧,\alpha }{}}\overline{\psi }_{𝐧,\alpha }(n_0)\psi _{𝐧,\alpha }(n_0)`$
$`\widehat{x}^+\widehat{x}`$ $`=`$ $`{\displaystyle \underset{𝐧,\alpha }{}}\widehat{x}_{𝐧,\alpha }^+\widehat{x}_{𝐧,\alpha },`$
$`\widehat{x}^+x(n_0)`$ $`=`$ $`{\displaystyle \underset{𝐧,\alpha }{}}\widehat{x}_{𝐧,\alpha }^+x_{𝐧,\alpha }(n_0)`$ (3)
and the quark action is

$`S_F`$ $`=`$ $`{\displaystyle \underset{m_0,n_0}{}}\overline{\psi }(m_0)Q(m_0,n_0)\psi (n_0)`$ (4)
$`=`$ $`{\displaystyle \underset{m,\alpha ,n,\beta }{}}\overline{\psi }_{𝐦,\alpha }(m_0)Q_{𝐦\alpha ,𝐧,\beta }(m_0,n_0)\psi _{𝐧,\beta }(n_0).`$
Let us come back to the evaluation of the quark determinant in a given baryonic sector.

The first step is to write the unconstrained determinant as a trace in Fock space

$$detQ=Tr\widehat{𝒯}.$$
(5)
$`\widehat{𝒯}`$ is the transfer matrix whose definition will be given in the next Section. The second step is to impose the restriction to a sector with baryon number $`n_B`$ by inserting in the trace the appropriate projection operator $`P_{n_B}`$

$$detQ|_{n_B}=Tr\left(\widehat{𝒯}P_{n_B}\right).$$
(6)
To rewrite the above expression as a determinant we choose coherent states  as a basis

$$|x>=|\mathrm{exp}(x\widehat{x}^+)>,|y>=|\mathrm{exp}(y\widehat{y}^+)>.$$
(7)
In such a basis the trace takes again the form of a Berezin integral

$`detQ|_{n_B}`$ $`=`$ $`{\displaystyle [dx^+dxdy^+dy]\mathrm{exp}\left(x^+xy^+y\right)<x,y|\widehat{𝒯}P_{n_B}|x,y>}`$
which will be expressed in terms of the determinant of a modified quark matrix.

So far for the strategy. Now we show how we get a series expansion, what is its the physical interpretation and what is the final result.

The kernel $`<x,y|\widehat{𝒯}P_{n_B}|x,y>`$ has the integral form

$`<x,y|\widehat{𝒯}P_{n_B}|x,y>`$ $`=`$ $`{\displaystyle }[dz^+dzdw^+dw]\mathrm{exp}(z^+zw^+w)`$ (9)
$`<x,y|\widehat{𝒯}|z,w><z,w|P_{n_B}|x,y>.`$
The kernel $`<x,y|\widehat{𝒯}|z,w>`$ will be given in the next Section, while that of $`P_{n_B}`$ is immediately calculated

$`<z,w|P_{n_B}|x,y>`$ $`=`$ $`{\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}(1)^{n_B}{\displaystyle \frac{1}{((n_B+r)!r!)^2}}`$ (10)
$`<(\widehat{y}w^+)^r(\widehat{x}z^+)^{(n_B+r)}(x\widehat{x}^+)^{n_B+r}(y\widehat{y}^+)^r>.`$
Since $`\widehat{x}^+,\widehat{y}^+`$ are creation operators of quarks and antiquarks respectively, we see that the r-th term of this series gives the gauge-invariant contribution of $`n_B`$ valence quarks plus $`r`$ quark-antiquark pairs.

Needless to say, for $`n_B=0`$, $`detQ|_{n_B}`$ does not reduce to the unconstrained determinant: Indeed also baryonic states are present in the unconstrained determinant, while they are not in $`detQ|_{n_B=0}`$. In QCD at nonvanishing temperature it makes a difference whether we impose or not the condition $`n_B=0`$. In view of the relatively low value of the critical temperature with respect to the nucleon mass, however, we do not expect significant effects from this restriction unless we go to exceedingly high tempreatures.

Let us now proceed to derive our final result. By evaluating the vacuum expectation values appearing in the last equation we express the kernel of the projection operator in terms of Grassmann variables only

$`<z,w|P_{n_B}|x,y>={\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}(1)^{n_B}{\displaystyle \frac{1}{(n_B+r)!r!}}(z^+x)^{n_B+r}(w^+y)^r.`$ (11)
To evaluate the integral of Eq.(9) we rewrite the above equation in exponential form

$`<z,w|P_{n_B}|x,y>={\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(n_B+r)!r!}}{\displaystyle \frac{^{n_B+r}}{j_1^{n_B+r}}}{\displaystyle \frac{^r}{j_2^r}}\mathrm{exp}(j_1z^+xj_2w^+y)|_{j_1=j_2=0}.`$ (12)
The integrals of Eqs.(2),(9) are gaussian and we get the constrained determinant in terms of the determinant of a modified quark matrix

$$detQ|_{n_B}=\underset{r=0}{\overset{\mathrm{}}{}}\frac{1}{(n_B+r)!r!}\frac{^{n_B+r}}{j_1^{n_B+r}}\frac{^r}{j_2^r}\left(Q+\delta Q\right)|_{j_1=j_2=0}.$$
(13)
The variation of the quark matrix is

$`\delta Q`$ $`=`$ $`K[(j_11)(1+\gamma _0)d(N_01)U_0(N_01)T_0^{(+)}+`$ (14)
$`(j_21)(1\gamma _0)d(N_0)T_0^{()}U_0^+(N_01)`$
$`+(j_1j_21)d(N_0)C(N_0)(1+\gamma _0)],`$
where

$`\left(T^{(\pm )}\right)_{m,n}`$ $`=`$ $`\delta _{m_0\pm 1,n_0}\delta _{𝐦,𝐧},`$
$`\left(d(N_0)\right)_{m,n}`$ $`=`$ $`\delta _{m_0,N_0}\delta _{m_0,n_0}\delta _{𝐦,𝐧},`$ (15)
$`K`$ is the hopping parameter and the matrix $`C`$ is given in the next Section. The particular time $`N_0`$ appearing depends on the fact that we inserted the projection operator at the latest time, but the value of $`detQ|_{n_B}`$ does not depend on this arbitrary choice, as it will be clear from the derivation.

## 3 Evaluation of the quark matrix in a given baryonic sector

In the first part of this Section we report the definition of the transfer matrix $`\widehat{𝒯}`$ and the proof of Eq. (5). These results have been derived by Lüscher , and we will follow his paper even in the notation with minor changes. Then we will evaluate the integral (9) getting the quark matrix in a given baryonic sector.

To start with we report the explicit expression of the quark matrix

$$Q=K\left[(1+\gamma _0)U_0T^{(+)}+(1\gamma _0)T^{()}U_0^++2C\right]B$$
(16)
where

$`C`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{3}{}}}\left(U_jU_j^+\right)\gamma _j,`$
$`B`$ $`=`$ $`11K{\displaystyle \underset{j=1}{\overset{3}{}}}\left(U_j+U_j^+\right)`$ (17)
and the lattice spacing has been set equal to 1. The modifications necessary at nonzero temperature are obvious. The transfer matrix is defined in terms of the operator

$`\widehat{T}_F(n_0)`$ $`=`$ $`\mathrm{exp}\left(2K\widehat{x}^+B(n_0)^{\frac{1}{2}}c(n_0)B(n_0)^{\frac{1}{2}}\widehat{y}^+\right)`$ (18)
$`\mathrm{exp}\left(\widehat{x}^+M(n_0)\widehat{x}\widehat{y}^+M(n_0)\widehat{y}\right)`$
according to the ordered product

$$\widehat{𝒯}=𝒥\underset{n_0=N_0}{\overset{N_01}{}}\left(\widehat{T}_F(n_0)\right)^+\widehat{T}_F(n_0+1).$$
(19)
$`𝒥`$ is the jacobian of a transformation which will be defined later, and two new matrices have made their appearance 111Lüscher assumes the gauge $`U_0=11`$ to prove reflection positivity. We are not concerned here with this most important issue and do not fix the gauge.

$`M`$ $`=`$ $`\mathrm{ln}\left((2K)^{\frac{1}{2}}U_0^+B^{\frac{1}{2}}\right)`$
$`c`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{3}{}}}\left(U_jU_j^+\right)i\sigma _j.`$ (20)
The $`\sigma _j`$ are the Pauli matrices and the operator $`\widehat{T}_F(n_0)`$ depends on the time $`n_0`$ only throu the dependence on it of the gauge fields.

We then use the following equations. Firstly, for arbitrary matrices $`M`$ and $`N`$

$$<x|\mathrm{exp}\widehat{x}^+M\widehat{x}|x^{}>=\mathrm{exp}\left(x^+\mathrm{exp}Mx^{}\right),$$
(21)
$$<x|\mathrm{exp}(\widehat{x}^+M\widehat{x})\mathrm{exp}(\widehat{x}^+N\widehat{x})|x^{}>=\mathrm{exp}\left(x^+e^Me^Nx^{}\right).$$
(22)
Secondly, if $`B=B(\widehat{x}^+)`$ and $`C=C(\widehat{x})`$ are operators which depend on $`\widehat{x}^+,\widehat{x}`$ only, for any operator A

$$<x|BAC|x^{}>=B(x^+)A(x^+,x^{})C(x^{}).$$
(23)
We thus get the kernel

$`<x(n_0),y(n_0)|\widehat{T}_F^+(n_0)\widehat{T}_F(n_0+1)|x(n_0+1),y(n_0+1)>=`$
$`\mathrm{exp}\left(2Kx(n_0)B(n_0)^{\frac{1}{2}}c(n_0)B(n_0)^{\frac{1}{2}}y(n_0)\right)`$
$`\mathrm{exp}\left(x^+(n_0)e^{M(n_0)}e^{M(n_0+1)}x(n_0+1)y^+(n_0)e^{M(n_0)}e^{M(n_0+1)}y(n_0+1)\right)`$
$`\mathrm{exp}\left(2Ky(n_0+1)B(n_0+1)^{\frac{1}{2}}c(n_0+1)B(n_0+1)^{\frac{1}{2}}x(n_0+1)\right).`$ (24)
Now to to go back from the Grassmann variables $`x`$’ and $`y`$’ to the quark field $`\psi `$ we only need the relations

$`\psi `$ $`=`$ $`B^{\frac{1}{2}}\left(\begin{array}{c}x\\ y^+\end{array}\right)`$
$`\overline{\psi }`$ $`=`$ $`\left(\begin{array}{cc}x^+& y\end{array}\right)\gamma _0B^{\frac{1}{2}}.`$ (25)
This transformation cancels out the jacobian $`𝒥`$. Inserting the identity between the factors in Eq.(19), using the equality

$$\overline{\psi }_n\left(B\psi \right)_n=x_n^+x_n+y_n^+y_n,$$
(26)
and collecting all the exponents we get Eq. (5). At this point two observations are in order. The first is that the interpretation of $`\widehat{x},\widehat{y}`$ as annihilation operators of quarks and antiquarks follows from the action of the charge conjugation on the quark field

$`\psi ^{}`$ $`=`$ $`𝒞^1\overline{\psi }`$
$`\overline{\psi }^{}`$ $`=`$ $`\stackrel{~}{𝒞}\overline{\psi },`$ (27)
where, with Lüscher’s convention for the $`\gamma `$-matrices,

$$𝒞=\gamma _0\gamma _2.$$
(28)
In fact these tranformations imply

$`x^{}`$ $`=`$ $`i\sigma _2y`$
$`y^{}`$ $`=`$ $`i\sigma _2x.`$ (29)
The second observation is that all the above construction requires the Wilson parameter $`r`$ to be equal to 1.

There remains the evaluation of Eq. (9). We notice that

$$<x,y|\widehat{𝒯}|z,w>=\tau \mathrm{exp}\left(x^+(N_01)Dz+y^+(N_01)D^{}w+wEz\right),$$
(30)
where $`\tau `$ does not depend on $`z,w`$, the star means complex conjugation and

$`D`$ $`=`$ $`B^{\frac{1}{2}}(N_01)U_0(N_01)B^{\frac{1}{2}}(N_0)`$
$`E`$ $`=`$ $`B^{\frac{1}{2}}(N_0)c(N_0)B^{\frac{1}{2}}(N_0).`$ (31)
We must therefore perform the gaussian integral

$`<x,y|\widehat{𝒯}P_{n_B}|x,y>={\displaystyle }[dz^+dzdw^+dw]\tau \mathrm{exp}(z^+zw^+w)`$
$`\mathrm{exp}\left(j_1z^+xj_2w^+y+x^+(N_01)Dz+y^+(N_01)D^{}w+wEz\right)`$
$`=\tau \mathrm{exp}\left(j_1x^+(N_01)Dxj_2y^+(N_01)D^{}y+j_1j_2yEx\right).`$ (32)
Comparing the kernel of the unconstrained transfer matrix with the above and taking into account the antiperiodic boundary conditions of the quark field in euclidean time we get the expression of Eq. (14) for $`\delta Q`$.

## 4 Conclusion

In principle by the present approach one can perform investigations of QCD in any definite baryonic sector at any temperature, starting from baryon number zero up to high baryonic densities. In practice, apart from the lowest baryonic numbers, we are confined to small physical volumes. With the volumes actually accessible in numerical simulations, perhaps this is not a too severe limitation, and we can hope to get reasonable results with the first few terms at least for not too high density and temperature.

One could think of different applications of our expansion by further specifying or changing the projection operator. In the case of baryonic number 1, for instance, by evaluating the expectation value of the total angular momentum  one can disentangle the valence and cloud contributions to the spin of the nucleon. Remaining to the lowest baryonic numbers, one might try to attack the study of the deuteron.

Another obvious field where our approach is potentially relevant is the QCD phase diagram. In particular one could investigate the various superconductive phases with spontaneous breaking of the color symmetry (see for instance ) by replacing the projection operator by the trial state of the BCS theory.

### Acknowledgements

It is a pleasure to thank dr.G.DiCarlo for many conversations about QCD at finite baryonic density.


text:
# Appendix A Appendix

## Appendix A Appendix

In this appendix we will provide some details on the derivation of the Killing spinor equations and conditions (14)-(22).

As already noted, the only consistent way to obtain a consistent set of differential equations from the supersymmetry variations (2) and (3) is to impose projections (12) on the spinors. We also provide for convenience the expressions for $`P_{\mu ij}`$ and $`Q_{\mu ij}`$ defined in (5). For the diagonal matrix $`L_\alpha ^i`$ in (8), it is convenient to represent them as forms in the index $`\mu `$,

$$P_{ij}=\left(\begin{array}{ccc}\lambda _1& gA^3\mathrm{sinh}\lambda _{12}& gA^2\mathrm{sinh}\lambda _{31}\\ gA^3\mathrm{sinh}\lambda _{12}& \lambda _2& gA^1\mathrm{sinh}\lambda _{23}\\ gA^2\mathrm{sinh}\lambda _{31}& gA^1\mathrm{sinh}\lambda _{23}& \lambda _3\end{array}\right)$$
(A.1)
and

$$Q_{ij}=\left(\begin{array}{ccc}0& gA^3\mathrm{cosh}\lambda _{12}& gA^2\mathrm{cosh}\lambda _{31}\\ gA^3\mathrm{cosh}\lambda _{12}& 0& gA^1\mathrm{cosh}\lambda _{23}\\ gA^2\mathrm{cosh}\lambda _{31}& gA^1\mathrm{cosh}\lambda _{23}& 0\end{array}\right).$$
(A.2)
We start with the $`i=1`$ case in the gaugino equation, $`\delta \chi _1=0`$, which implies two different equations: factorizing $`\widehat{\mathrm{\Gamma }}_2\widehat{\mathrm{\Gamma }}_3ϵ`$ we get

$$\frac{d\lambda _1}{d\rho }+\frac{2}{3}\frac{d\mathrm{\Phi }}{d\rho }=e^{\mathrm{\Phi }+\lambda _1}\frac{F_{23}^1}{\alpha _2\alpha _3}\frac{g}{4}e^\mathrm{\Phi }(e^{2\lambda _1}e^{2\lambda _2}e^{2\lambda _3}),$$
(A.3)
where $`F_{23}^1`$ is defined in (19). In addition, from terms proportional to $`\widehat{\mathrm{\Gamma }}_2\mathrm{\Gamma }_3ϵ`$ we get

$$e^{\mathrm{\Phi }+\lambda _1}\frac{F_{\rho 1}^1}{\alpha _1}+g\left(\frac{A_3^3}{\alpha _3}\mathrm{sinh}\lambda _{12}\frac{A_2^2}{\alpha _2}\mathrm{sinh}\lambda _{31}\right)=0,$$
(A.4)
where $`F_{\rho 1}^1=_\rho A_1^1`$, or equivalently, after the change of variables (29) (and setting $`g=2`$)

$$\frac{b_1\dot{c}_1}{a_1}+\frac{c_3}{a_3}\left(\frac{b_1}{b_2}\frac{b_2}{b_1}\right)\frac{c_2}{a_2}\left(\frac{b_3}{b_1}\frac{b_1}{b_3}\right)=0.$$
(A.5)
The four additional equations corresponding to $`\delta \chi _i=0`$, for $`i=2,3`$ can be obtained from (A.3) and (A.4) by cyclic permutations in the indices $`1,2,3`$. Using then the constraint $`\lambda _1+\lambda _2+\lambda _3=0`$ we get

$$\frac{d\lambda _1}{d\rho }=\frac{e^\mathrm{\Phi }}{3}\left(2\frac{e^{\lambda _1}}{\alpha _2\alpha _3}F_{23}^1\frac{e^{\lambda _2}}{\alpha _3\alpha _1}F_{31}^2\frac{e^{\lambda _3}}{\alpha _1\alpha _2}F_{12}^3\right)\frac{g}{6}e^\mathrm{\Phi }(2e^{2\lambda _1}e^{2\lambda _2}e^{2\lambda _3})$$
(A.6)
and

$$\frac{d\mathrm{\Phi }}{d\rho }=\frac{1}{2}e^\mathrm{\Phi }\left(\frac{e^{\lambda _1}}{\alpha _2\alpha _3}F_{23}^1+\frac{e^{\lambda _2}}{\alpha _3\alpha _1}F_{31}^2+\frac{e^{\lambda _3}}{\alpha _1\alpha _2}F_{12}^3\right)+\frac{g}{8}e^\mathrm{\Phi }(e^{2\lambda _1}+e^{2\lambda _2}+e^{2\lambda _3}).$$
(A.7)
Next we turn to the gravitino variation, $`\delta \psi _\mu =0`$, where we should distinguish two possibilities, according to whether $`\mu `$ is a coordinate on the wrapped 3-sphere or on the unwrapped directions. Thus, if $`\mu =\sigma _1`$, from $`\widehat{\mathrm{\Gamma }}_2\widehat{\mathrm{\Gamma }}_3ϵ`$ we get

$$\omega _1^{23}+gA_1^1\mathrm{cosh}\lambda _{23}\alpha _1\frac{e^\mathrm{\Phi }}{6}\left(\frac{e^{\lambda _2}}{\alpha _2}F_{\rho 2}^2+\frac{e^{\lambda _3}}{\alpha _3}F_{\rho 3}^35\frac{e^{\lambda _1}}{\alpha _1}F_{\rho 1}^1\right)=0$$
(A.8)
and, from terms proportional to $`\widehat{\mathrm{\Gamma }}_2\mathrm{\Gamma }_3ϵ`$

$$\frac{1}{\alpha _1}\frac{d\alpha _1}{d\rho }\frac{1}{6}e^\mathrm{\Phi }\left(\frac{e^{\lambda _1}}{\alpha _2\alpha _3}F_{23}^15\frac{e^{\lambda _2}}{\alpha _3\alpha _1}F_{31}^25\frac{e^{\lambda _3}}{\alpha _1\alpha _2}F_{12}^3\right)\frac{g}{24}e^\mathrm{\Phi }(e^{2\lambda _1}+e^{2\lambda _2}+e^{2\lambda _3})=0,$$
(A.9)
where we have used the spin connection component $`\omega _i^{i\rho }=\frac{d\alpha _i}{d\rho }`$ for the metric ansatz (9). The equations obtained from cyclicity of (A.8) and (A.9) correspond to the choices $`\mu =\sigma _2,\sigma _3`$. The generalized twist (14) is then derived from (A.4) and (A.8). It amounts to turning on a gauge field given by (15).

It is important to verify that substituting back (15) into (A.4) and after using (A.6), (A.7) and (A.9), gives no new constraints for the various functions. After a straightforward but lengthly computation, one can show that this is indeed the case.

If $`\mu `$ is a coordinate on the unwrapped part of the worldvolume, using $`\omega _\mu ^{\mu \rho }=e^f\frac{df}{d\rho }`$ and comparing terms proportional to $`\mathrm{\Gamma }_\mu \mathrm{\Gamma }_7ϵ`$ we get and equation relating the warp factor and the dilaton as

$$\frac{df}{d\rho }=\frac{1}{3}\frac{d\mathrm{\Phi }}{d\rho }.$$
(A.10)
This then leads to (20) in the main text (a possible constant of integration can be absorbed into a rescaling of the corresponding unwrapped coordinates). Also from comparison of terms proportional to $`\mathrm{\Gamma }_\mu \mathrm{\Gamma }_7\mathrm{\Gamma }_1\widehat{\mathrm{\Gamma }}_1ϵ`$ we obtain a constraint for the field strength

$$e^{\lambda _1}\frac{F_{\rho 1}^1}{\alpha _1}+e^{\lambda _2}\frac{F_{\rho 2}^2}{\alpha _2}+e^{\lambda _3}\frac{F_{\rho 3}^3}{\alpha _3}=0,$$
(A.11)
which holds identically from (A.4). Finally, considering the gravitino equation for $`x^\mu =\rho `$ gives after some algebra a simple differential equation for the spinor $`ϵ`$ with solution given by (22).


text:
# The efficient generation of unstructured control volumes in 2D and 3D

## 1 Introduction

An unstructured mesh generator can solve many of the problems associated with the structured meshes. The unstructured meshes are suitable for complex geometries, especially when we use the control volume method . The control volume method on the one hand, is based upon the simple physical principle of a flux balance. Nevertheless, the generated meshes have to fulfill some restrictions for such method. One of the conditions says than lines connecting the central point in the control volume with neighbour central points have to be perpendicular to the sides of the neighbour volumes. However, the neighbourhood of each neighbour or cell in an unstructured mesh must be defined explicitly. This is one disadvantage of unstructured meshes because it has to reserve large storage of the computer memory. Nevertheless, the advantages like easy for handling adaptability in time, ability to generate meshes about arbitrary geometries. We can divide such meshes into the cell-vertex control volumes and the cell-centered control volumes. In the cell-vertex volumes, all elements containing the relevant point are applied as a control volume, so the control volumes overlap. In the cell-centered control volumes the elements are subdivided and the control volumes do not overlap.

In this paper we propose a generator for unstructured cell-centered volumes, which is useful for the control volume method . The mesh is constructed for a number of points randomly located in a domain and the domain border.

## 2 Mathematical background

We turn our attention in two and three dimensional space. Let us consider a domain $`\mathrm{\Omega }R^3`$ with smooth its boundary $`\mathrm{\Omega }`$. In such domain and boundary we generate randomly a set of points $`P_i(a,b,c)`$, where $`i=1..N`$. The points establish cell-centers, around of which we determine convex polygons. The polygons are the non-overlapping control volumes which are applied in the control volume method . We can formulate two groups of basic assumptions necessary for the meshes construction as follows:

1. global conditions:

such domain has to be consistent,

edges of some polygon have the same lengths to the corresponding edges of neighbour polygons,
2. local conditions:

inside a polygon one may found only one random point called the cell-center or the point,

an edge of the polygon has to be perpendicular to the line connecting two points.

If we fulfill the global assumptions, we may discretize the domain $`\mathrm{\Omega }`$ through its division into the triangles in $`R^2`$ space or into the tetrahedrons in $`R^3`$ space. Such discretization is a standard solution, which one may find in convex geometry  through Delaunay triangulation . Regarding to the global and local conditions we take into consideration a fact, that we cannot generate polygons in $`R^2`$ which overlap neighbour polygons. The fact also exist for the generated tetrahedrons in $`R^3`$ space.

In this paper we propose a novel algorithm for generation of control volumes. We start the algorithm in two dimensions because one can see how it works. We consider circles with unknown radiuses which centers are located in the points, which are randomly located in the domain $`\mathrm{\Omega }`$. Let us assume, that within the triangle created by the points the circles only in the one point crosses. The point is a vertice of control volumes generated around the point. Following that we need to solve the system of equations

$$(\widehat{x}a_l)^2+(\widehat{y}b_l)^2=r_l^2,(l=1,2,3).$$
(1)
Assuming, that the radiuses $`r_l`$ of the circles are unknown, we obtain a point $`Q(\widehat{x}(r_1,r_2,r_3),\widehat{y}(r_1,r_2,r_3))`$ which coordinates are solution of the system (1)

$$\begin{array}{c}\widehat{x}(r_1,r_2,r_3)=\frac{(b_2b_3)r_1^2+(b_1b_3)r_2^2(b_1b_2)r_3^2+(a_1^2+b_1^2)(b_2b_3)}{2\left[(a_1a_2)(b_1b_3)(a_1a_3)(b_1b_2)\right]}+\\ \frac{(a_2^2+b_2^2)(b_1b_3)+(a_3^2+b_3^2)(b_1b_2)}{2\left[(a_1a_2)(b_1b_3)(a_1a_3)(b_1b_2)\right]}\\ \widehat{y}(r_1,r_2,r_3)=\frac{(a_2a_3)r_1^2(a_1a_3)r_2^2+(a_1a_2)r_3^2(a_1^2+b_1^2)(a_2a_3)}{2\left[(a_1a_2)(b_1b_3)(a_1a_3)(b_1b_2)\right]}+\\ +\frac{(a_2^2+b_2^2)(a_1a_3)(a_3^2+b_3^2)(a_1a_2)}{2\left[(a_1a_2)(b_1b_3)(a_1a_3)(b_1b_2)\right]}\end{array}$$
(2)
The indexes in formulas (1) and (2) correspond to the local case.

In the global case we have random points $`P_i(a,b,c),i=1..N`$ which are generated inside the domain $`\mathrm{\Omega }`$ and its boundary $`\mathrm{\Omega }`$. We use Delaunay triangulation  that to find neighbourhood of point $`P_i`$ defined by several points $`P_{j(i,1)}\mathrm{}P_{j(i,k)}\mathrm{}P_{j(i,M(i))}`$. The index $`j(i,k)`$ is a function which establishes a point number having neighbourhood with the $`i`$-th point. The temporal index $`k`$ varies from $`1`$ to $`M(i)`$. The $`M(i)`$ is a number of points corresponding to the $`i`$-th point. Three points $`P_i,P_{j(i,k)},P_{j(i,k+1)}`$ define a triangle. When the temporal index $`k+1`$ exceeds $`M(i)`$ then we have $`P_{j(i,k+1)}=P_{j(i,1)}`$. Moreover, we take into consideration a condition $`j(i,k)i`$. We extend a definition of formula (2) for the global case putting $`(a_i,b_i)`$, $`(a_{j(i,k)},b_{j(i,k)})`$, $`(a_{j(i,k+1)},b_{j(i,k+1)})`$ and $`(r_i,r_{j(i,k)},r_{j(i,k+1)})`$ respectively into variables with established indexes. In our approach, we are looking for the radiuses $`r_i`$. Taking into consideration of such fact we formulate a range of radiuses variation as

$$\underset{k}{\mathrm{max}}\left[(L_{i,j(i,k)}r_{maxj(i,k)})\right]<r_i<r_{maxi},$$
(3)
where $`L_{i,j(i,k)}`$ is a distance between two neighbour points defined as

$$L_{i,j(i)}=\sqrt{(a_ia_{j(i)})^2+(b_ib_{j(i)})^2}.$$
(4)
Maximal radius of a circle is formed as

$$r_{maxi}=\underset{k}{\mathrm{min}}\left[h_{i,j(i,k)}\right],$$
(5)
where $`h_{i,j(i,k)}`$ is a distance defined in acute triangles as

$$\begin{array}{c}h_{i,j(i,k)}=\frac{|(b_{j(i,k)}b_{j(i,k+1)})a_i(a_{j(i,k)}a_{j(i,k+1)})b_i|}{\sqrt{(b_{j(i,k)}b_{j(i,k+1)})^2+(a_{j(i,k)}a_{j(i,k+1)})^2}}+\\ +\frac{|(a_{j(i,k)}a_{j(i,k+1)})b_{j(i,k+1)}(b_{j(i,k)}b_{j(i,k+1)})a_{j(i,k+1)}|}{\sqrt{(b_{j(i,k)}b_{j(i,k+1)})^2+(a_{j(i,k)}a_{j(i,k+1)})^2}}\end{array}.$$
(6)
In acute angles of right and obtuse triangles we have

$$h_{i,j(i,k)}=\underset{k}{\mathrm{min}}[L_{i,j(i,k)},L_{i,j(i,k+1)}].$$
(7)
The symbol $`r_{maxj(i,k)}`$ represents a maximal radius of the $`j`$-th point being a neighbour of the $`i`$-th point.

Using extended formula (2) and additional expressions (3)$`÷`$(7) one can construct an unstructured mesh. Fig. 1 shows an idea of the control volume generation using the expression (2) and the limit (3). We can observe that the local condition - an edge of the polygon has to be perpendicular to the line connecting two points - is fulfilled. The general algorithm uses some optimisation that to find the points (2) being the vertices of the polygons (control volumes). In this paper we neglect the construction of an objective function because it is not the subject of our consideration. We applied Rosenbrock’s method  for the local solution. For the global solution we use the soft selection method  which is a global optimisation method. Prediction of the unstructured mesh under such way is not only the one way of the mesh construction. It is possible to construct another type of mesh without restriction that expression (2) is a solution of eqn. (1). But we have to restrict the expression (3) because it guaranties right mesh construction for the control volume method . Let us assume that a radius of the circle overlaps other radiuses of the neighbour circles. Of course we have to restrict the condition (3). Corresponding to previous and present mesh construction we have

$$L_{i,j(i,k)}\left(r_i+r_{j(i,k)}\right)0.$$
(8)
Following eqn. (8) the point coordinates given by expression (2) are not the solution of the system (1) but they become the coordinates of the point predicted by lines connecting two points of crossing circles. From the other hand, we can assume that a radius of the circle does not overlap other radiuses of the neighbour circles. In such situation we have

$$L_{i,j(i,k)}\left(r_i+r_{j(i,k)}\right)>0.$$
(9)
One can say, that conditions (8) and (9) allow us to obtain different forms of unstructured meshes. However, in local solution is possible some combination of expressions (8) and (9) respectively.

Similar to previous results we consider a set of random points $`P_i(a,b,c)`$ in $`R^3`$ space. Analogously to the previous considerations, we use spatial Delaunay triangulation that to find neighbourhood of point $`P_i`$ which is defined by several points $`P_{j(i,k,l)}`$. The index $`j(i,k,l)`$ is a function which establishes a point number having neighbourhood with the $`i`$-th point. The temporal index $`k`$ describes a number of tetrahedrons having the common $`i`$-th point. It varies from $`1`$ to $`N(i)`$. The next index $`l`$ varies from $`1`$ to $`3`$ and it is a point number corresponding to the $`k`$-th tetrahedron. The four points $`P_i,P_{j(i,k,1)},P_{j(i,k,2)},P_{j(i,k,3)}`$ define a $`k`$-th tetrahedron. We search radiuses of spheres which centers are located in the points. In each tetrahedron crosses the four radiuses and they create a point which is the vertice of the control volume. Solving the following system

$$(\widehat{x}a_l)^2+(\widehat{y}b_l)^2+(\widehat{z}c_l)^2=r_l^2,(l=\mathrm{1..4}),$$
(10)
we can find the local vertices coordinates

$$\widehat{x}(r_1,r_2,r_3,r_4)=\frac{W_x}{W},\widehat{y}(r_1,r_2,r_3,r_4)=\frac{W_y}{W},\widehat{z}(r_1,r_2,r_3,r_4)=\frac{W_z}{W},$$
(11)
where

$`W=\left|\begin{array}{ccc}2(a_1a_2)& 2(b_1b_2)& 2(c_1c_2)\\ 2(a_1a_3)& 2(b_1b_3)& 2(c_1c_3)\\ 2(a_1a_4)& 2(b_1b_4)& 2(c_1c_4)\end{array}\right|,W_x=\left|\begin{array}{ccc}\mathrm{}_1& 2(b_1b_2)& 2(c_1c_2)\\ \mathrm{}_2& 2(b_1b_3)& 2(c_1c_3)\\ \mathrm{}_3& 2(b_1b_4)& 2(c_1c_4)\end{array}\right|,`$
$`W_y=\left|\begin{array}{ccc}2(a_1a_2)& \mathrm{}_1& 2(c_1c_2)\\ 2(a_1a_3)& \mathrm{}_2& 2(c_1c_3)\\ 2(a_1a_4)& \mathrm{}_3& 2(c_1c_4)\end{array}\right|,W_z=\left|\begin{array}{ccc}2(a_1a_2)& 2(b_1b_2)& \mathrm{}_1\\ 2(a_1a_3)& 2(b_1b_3)& \mathrm{}_2\\ 2(a_1a_4)& 2(b_1b_4)& \mathrm{}_3\end{array}\right|`$
and

$`\mathrm{}_1=r_2^2r_1^2+a_1^2a_2^2+b_1^2b_2^2+c_1^2c_2^2,`$
$`\mathrm{}_2=r_3^2r_1^2+a_1^2a_3^2+b_1^2b_3^2+c_1^2c_3^2,`$
$`\mathrm{}_3=r_4^2r_1^2+a_1^2a_4^2+b_1^2b_4^2+c_1^2c_4^2.`$
Similar to previous mesh construction, we define lower and upper limits of the radiuses variation

$$\underset{k,l}{\mathrm{max}}\left[(h_{i,j(i,k,l)}r_{maxj(i,k)})\right]<r_i<r_{maxi}.$$
(14)
The symbol $`h_{i,j(i,k,l)}`$ represents height of a triangle being the wall of the tetrahedron and for the acute triangles is defined by the following formula

$$\begin{array}{c}h_{i,j(i,k,l)}=\sqrt{\frac{\left(BB_{12}CC_{32}BB_{32}CC_{12}\right)^2\left(AA_{32}CC_{12}AA_{12}CC_{32}\right)^2+}{AA_{32}^2+BB_{32}^2+CC_{32}^2}}+\\ \sqrt{\frac{\left(AA_{12}BB_{32}AA_{32}BB_{12}\right)^2}{AA_{32}^2+BB_{32}^2+CC_{32}^2}}\end{array},$$
(15)
where

$`AA_{12}`$ $`=`$ $`a_ia_{j(i,k,l)},AA_{32}=a_{j(i,k,l+1)}a_{j(i,k,l)},`$
$`BB_{12}`$ $`=`$ $`b_ib_{j(i,k,l)},BB_{32}=b_{j(i,k,l+1)}b_{j(i,k,l)},`$
$`CC_{12}`$ $`=`$ $`c_ic_{j(i,k,l)},CC_{32}=c_{j(i,k,l+1)}c_{j(i,k,l)}.`$
In acute angles of right and obtuse spatial triangles we have

$$h_{i,j(i,k,l)}=\underset{l}{\mathrm{min}}[L_{i,j(i,k,l)},L_{i,j(i,k,l+1)}],$$
(16)
where $`L_{i,j(i,k,l)}`$ is a spatial distance defined between two points as

$$L_{i,j(i,k,l)}=\sqrt{\left(a_ia_{j(i,k,l)}\right)^2+\left(b_ib_{j(i,k,l)}\right)^2+\left(c_ic_{j(i,k,l)}\right)^2}.$$
(17)
When the temporal index $`l+1`$ exceeds $`3`$ then we put $`P_{j(i,k,l+1)}=P_{j(i,k,1)}`$. Maximal radius of a sphere is defined as

$$r_{maxi}=\underset{k}{\mathrm{min}}\left[H_{i,j(i,k)}\right].$$
(18)
However, the function $`H_{i,j(i,k)}`$ is height of a tetrahedron created by the vertices $`P_i,P_{j(i,k,1)},P_{j(i,k,2)},P_{j(i,k,3)}`$ and for the acute tetrahedron we have

$$H_{i,j(i,k)}=\frac{|Aa_i+Bb_i+Cc_i+D|}{\sqrt{A^2+B^2+C^2}}$$
(19)
where

$`A=m_{21}n_{31}m_{31}n_{21},`$ $`B=d_{31}n_{21}d_{21}n_{31},`$
$`C=d_{21}m_{31}d_{31}m_{21},`$ $`D=Aa_{j(i,k,1)}Bb_{j(i,k,1)}Cc_{j(i,k,1)}`$
and

$`d_{21}=a_{j(i,k,2)}a_{j(i,k,1)},m_{21}=b_{j(i,k,2)}b_{j(i,k,1)},n_{21}=c_{j(i,k,2)}c_{j(i,k,1)},`$
$`d_{31}=a_{j(i,k,3)}a_{j(i,k,1)},m_{31}=b_{j(i,k,3)}b_{j(i,k,1)},n_{31}=c_{j(i,k,3)}c_{j(i,k,1)}.`$
But for the right and obtuse tetrahedrons we have

$$H_{i,j(i,k)}=\underset{l}{\mathrm{min}}\left[h_{i,j(i,k,l)}\right].$$
(20)
The symbol $`r_{maxj(i,k)}`$ represents a maximal radius of $`j`$-th point being a neighbour of the $`i`$-th point.

## 3 Example of calculations

In practical mesh construction when expression (3) is fulfilled we obtain a unique solution. Fig. 2 shows an example of mesh created under restriction that formula (2) is a solution of the system (1). Fig. 2a presents triangulation of some points in which generated circles. Within the triangle crosses the circles only in one point. According to the previous description we draw control volumes presented by Fig. 2b. Nevertheless, we can construct a mesh without restriction that formula (2) is a solution of the system (1). We can apply conditions (8) or (9) respectively. In this case we generate a mesh much more easier than in previous case. Expressions (8) and (9) extend our considerations for different generation of meshes.

## 4 Concluding remarks

In this paper we elaborate a new method of generation of unstructured meshes. The algorithm is prepared for the control volume method in two- and three dimensional space. For random location of points we use Delaunay triangulation before the mesh generation. Our method bases in 3D on circles generation in the points. Within one triangle created by three points crosses the circles in one point. The point is a vertice of a control volume. We extend our considerations for the triangles in which are overlapping or non-overlapping circles. The extension allow us to generate meshes much more easier than in previous case. In opposite to the Thyssen polygons we can generate control volumes for the right and obtuse triangles. However, our method is suitable and easy to use in generation of control volumes in three dimensional space.