text:
# Density field in extended Lagrangian perturbation theory

## I Introduction

The Lagrangian approximation for structure formation in cosmological fluids provides a relatively accurate model even in quasi-linear regime, where density fluctuation becomes unity. The Zel’dovich approximation (hereafter ZA) zel ; buchert92 ; coles ; saco , a linear Lagrangian approximation for dust fluid, describes the evolution of density fluctuation better than Eulerian approximation munshi ; sahsha . Although ZA gives an accurate description until quasi-linear regime develops, ZA cannot describe the model after the formation of caustics. In ZA, even after the formation of caustics, the fluid elements keeps moving in the direction set up by the initial condition. ZA cannot describe compact and high density structures such as pancakes, skeletons, or clumps, while N-body simulation shows the presence of clumps with a very wide range in mass at any given time davis . In general, after the formation of caustics, hydrodynamical description becomes invalid.

In order to proceed with a hydrodynamical description without the formation of caustics, the qualitative pressure gradient zel82  and thermal velocity scatter kotok ; shandarin  in a collisionless medium had been discussed. After that, the ‘adhesion approximation’ gurbatov  was proposed from the consideration of nonlinear wave equations like Burgers’ equation. In the ‘adhesion approximation’, an artificial viscosity term is added to ZA. As another modification, a Gaussian cutoff is applied to the initial power spectrum of density fluctuation and then evolves with the ZA. This modified approximation is called the ‘Truncated Zel’dovich approximation’ (hereafter TZA) cms ; mps . These modified approximations successfully avoid the formation of caustics and describe the evolution for long time. However, the physical origin of the modification is not clarified.

We reconsider the basic, fundamental equation for the motion of matter. The collisionless Boltzmann equation BT  describes the motion of matter in phase space. The basic equations of hydrodynamics are obtained by integrating the collisionless Boltzmann equation over velocity space. In past approximations, such as ZA and its modified version, velocity dispersion was ignored. Buchert and Domínguez budo  argued that the effect of velocity dispersion become important beyond the caustics. They also argued that models for large-scale structure should rather be constructed for a flow which describes the average motion of a multi-stream system. Then they showed that when the velocity dispersion is still small and can be considered isotropic, that gives effective ’pressure’ or viscosity terms. Furthermore, they argued the relation between mass density $`\rho `$ and pressure $`P`$, i.e. an ‘equation of state’. If the relation between the density of matter and pressure seems barotropic, the equation of state should take the form $`P\rho ^{5/3}`$. Buchert et al. bdp  showed how the viscosity term is generated by the effective pressure of a fluid under the assumption that the peculiar acceleration is parallel to the peculiar velocity; Domínguez domi00 ; domi0106  clarified that a hydrodynamic formulation is obtained via a spatial coarse graining in a many-body gravitating system, and the viscosity term in the ‘adhesion approximation’ can be derived by the expansion of coarse-grained equations.

Domínguez domi0103  also reported on a study of the spatially coarse-grained velocity dispersion in cosmological N-body simulations. The analysis showed that polytrope exponent becomes $`\gamma 5/3`$ in a quasi-nonlinear regime, and $`\gamma 2`$ in a strongly nonlinear regime. Domínguez and Melott dome0310  discussed polytrope exponents of velocity dispersion in N-body simulations. According to their results, the exponents depend on a model of initial density fluctuation.

From these points, the extension of Lagrangian perturbation theory to cosmological fluids with pressure has been considered. Actually, Adler and Buchert adler  have formulated the Lagrangian perturbation theory for a barotropic fluid. Morita and Tatekawa moritate  and Tatekawa et al. tate02  solved the Lagrangian perturbation equations for a polytropic fluid up to second order for cases where the equations are solved easily. Hereafter, we call this model the ‘pressure model’.

In this paper, we analyze the density field which is described by the Lagrangian approximations; ZA, TZA, and first-order pressure model solutions. We calculate the cross-correlation function of density field between the Lagrangian approximation and N-body simulation. Furthermore we analyze the probability distribution function of density fluctuation for confirmation. From these results, we determine a polytrope exponent in the equation of state. From our analyses of the cross-correlation coefficient and probability distribution function of density fluctuation, we find that the value of the polytrope exponent seems to be $`5/3`$ for quasi-nonlinear evolution, as Buchert and Domínguez argued budo . However, for the determination of a proportion fixed number in equation of state, we must consider further physical processes or carry out high-resolutional N-body simulation.

This paper is organized as follows. In Sec. II, we present Lagrangian perturbative solutions: In Sec. II.1, we show the first-order solution of pressure model in the Einstein-de Sitter background. For comparison, in Sec. II.2 and II.3, we show the solution of ZA and the procedure of TZA.

In Sec. III, we compare the density field between the Lagrangian approximations and N-body simulation. In Sec.III.1, we calculate the cross-correlation coefficient of the density field. Though it seems that we can reach a conclusion in this analysis, it is insufficient. Therefore in Sec. III.2, we analyze the probability distribution function of density fluctuation. In Sec. IV, we discuss our results and state conclusions.

## II Lagrangian approximations in gravitational instability theory

### II.1 First-order solutions of the pressure model

In this section, we present perturbative solutions in the Lagrangian description. The matter model we consider is a self-gravitating fluid with mass density $`\rho `$ and ‘pressure’ $`P`$, which is given by the presence of velocity dispersion. The ‘pressure’ we adopt here is the same as was introduced by Buchert and Domínguez budo , i.e. the diagonal component of the velocity dispersion tensor when the velocity dispersion is assumed to be isotropic in the Jeans equation BT . In Lagrangian hydrodynamics, the comoving coordinates $`𝒙`$ of the fluid elements are represented in terms of Lagrangian coordinates $`𝒒`$ as

$$𝒙=𝒒+𝒔(𝒒,t),$$
(1)
where $`𝒒`$ are defined as initial values of $`𝒙`$, and $`𝒔`$ denotes the Lagrangian displacement vector due to the presence of inhomogeneities. From the Jacobian of the coordinate transformation from $`𝒙`$ to $`𝒒`$, $`Jdet(x_i/q_j)=det(\delta _{ij}+s_i/q_j)`$, the mass density is given exactly as

$$\rho =\rho _\mathrm{b}J^1.$$
(2)
We decompose $`𝒔`$ into the longitudinal and the transverse modes as $`𝒔=_𝒒S+𝑺^\mathrm{T}`$ with $`_𝒒𝑺^\mathrm{T}=0`$. In this paper, we show an explicit form of perturbative solutions in only Einstein-de Sitter universe. For a generic background universe, we obtained the perturbation solutions in our previous paper tate02 .

The transverse modes can be solved easily. The first-order solutions become as follows:

$$𝑺^\mathrm{T}a^0,a^{1/2}.$$
(3)
Because the solutions do not depend on ’pressure’, the solutions of this mode in ZA become the same in form. The transverse modes do not have a growing solution in a first-order approximation.

For the longitudinal modes, we carry out Fourier transformation with respect to the Lagrangian coordinates $`𝒒`$. If we assume a polytropic equation of state $`P=\kappa \rho ^\gamma `$ with a constant $`\kappa `$ and a polytrope exponent $`\gamma `$, we can write an explicit form of first-order perturbative solutions. In the Einstein-de Sitter background, the solutions are written in a relatively simple form. They are, for $`\gamma 4/3`$,

$$\widehat{S}(𝑲,a)a^{1/4}𝒥_{\pm 5/(86\gamma )}\left(\sqrt{\frac{2C_2}{C_1}}\frac{|𝑲|}{|43\gamma |}a^{(43\gamma )/2}\right),$$
(4)
where $`𝒥_\nu `$ denotes the Bessel function of order $`\nu `$, and for $`\gamma =4/3`$,

$$\widehat{S}(𝑲,a)a^{1/4\pm \sqrt{25/16C_2|𝑲|^2/2C_1}},$$
(5)
where $`C_14\pi G\rho _\mathrm{b}(a_{\mathrm{in}})a_{\mathrm{in}}^{3}/3`$ and $`C_2\kappa \gamma \rho _\mathrm{b}(a_{\mathrm{in}})^{\gamma 1}a_{\mathrm{in}}^{3(\gamma 1)}`$. $`\rho _b`$ and $`𝑲`$ mean background mass density and Lagrangian wavenumber, respectively. $`a_{\mathrm{in}}`$ means scale factor when an initial condition was given.

Here we notice the relation between the behaviors of the above solutions and the Jeans wavenumber, which is defined as

$$K_\mathrm{J}\left(\frac{4\pi G\rho _\mathrm{b}a^2}{\mathrm{d}P/\mathrm{d}\rho (\rho _\mathrm{b})}\right)^{1/2}.$$

The Jeans wavenumber, which gives a criterion for whether a density perturbation with a wavenumber will grow or decay with oscillation, depends on time in general. If the polytropic equation of state $`P=\kappa \rho ^\gamma `$ is assumed,

$$K_\mathrm{J}=\sqrt{\frac{3C_1}{C_2}}a^{(3\gamma 4)/2}.$$
(6)
Equation (6) implies that, if $`\gamma <4/3`$, $`K_\mathrm{J}`$ will be decreased and density perturbations with any scale will decay and disappear in evolution, and if $`\gamma >4/3`$, all density perturbations will grow to collapse. We rewrite the first-order solution Equation (4) with the Jeans wavenumber:

$$\widehat{S}(𝑲,a)a^{1/4}𝒥_{\pm 5/(86\gamma )}\left(\frac{\sqrt{6}}{|43\gamma |}\frac{|𝑲|}{K_\mathrm{J}}\right).$$
(7)
### II.2 Zel’dovich approximation

ZA was obtained as first-order solutions with dust fluid in Lagrangian description zel . The solutions are obtained from the solutions of the pressure model as a limit of weak pressure. For example, in E-dS model, when we take a limit $`\kappa 0`$ in Equations (4) and (5), the solutions converge to those of ZA:

$$\widehat{S}(𝑲,a)a,a^{3/2}.$$
(8)
ZA is known as perturbative solutions which describe the structure well in quasi-linear scale. However if the caustics appear, the solutions no longer have physical meaning.

### II.3 Truncated Zel’dovich approximation

During evolution, the small scale structure contracts and forms caustics. Therefore if we introduce some cutoff in the small scale, we will be able to avoid the formation of caustics cms ; mps . In TZA, for the avoidance of caustics, we introduce a Gaussian cutoff to the initial density spectrum as follows:

$$𝒫(k,t_{\mathrm{in}})𝒫(k,t_{\mathrm{in}})\mathrm{exp}\left(k^2/k_{NL}\right),$$
(9)
where $`k_{NL}`$ means ‘nonlinear wavenumber’, defined by

$$1=a(t)^2_{k_0}^{k_{NL}}𝒫(k,t_{\mathrm{in}})𝑑k.$$
(10)
The ‘nonlinear wavenumber’ depends on the scale factor. The relation between the Jeans wavenumber $`K_J`$ and the nonlinear wavenumber $`k_{NL}`$ will be discussed in Sec. IV.

## III Comparison between N-body simulation and Lagrangian approximations

In this section, we show a comparison between N-body simulation and Lagrangian approximations with two statistical quantities. In our previous paper tate02 , we showed that the effect of second order perturbation was still small just before shell-crossing. Therefore we consider only first order perturbation.

We analyze ZA zel , TZA cms ; mps , and the pressure model adler ; moritate ; tate02 . We establish the value of scale factor at $`z=0`$ with $`a=1`$. For the initial condition, we set the Gaussian density field with scale-free spectrum;

$$𝒫(k)k^n(n=1,0,1),$$
(11)
SCDM, and LCDM model. The initial condition was produced by COSMICS COSMICS .

For N-body simulation, we execute $`P^3M`$ code. The parameters of simulation were given as follows:

Number of particles $`:`$ $`N=64^3,N=128^3\text{(Figure }\text{3}\text{ only)},`$
Box size $`:`$ $`L=64h^1\text{Mpc},`$
Softening Length $`:`$ $`\epsilon =0.05h^1\text{Mpc},`$
Course-Graining Length $`:`$ $`l=1,2,4h^1\text{Mpc},`$
Hubble patameter $`:`$ $`h=0.71.`$
For CDM models, we choose cosmological parameters as follows:

SCDM $`:`$ $`\mathrm{\Omega }_m=1.0,\mathrm{\Omega }_\mathrm{\Lambda }=0.0,\sigma _8=0.84`$
LCDM $`:`$ $`\mathrm{\Omega }_m=0.27,\mathrm{\Omega }_\mathrm{\Lambda }=0.73,\sigma _8=0.84`$
In the pressure model, we choose a polytrope exponent $`\gamma =4/3,5/3`$: In the case of $`\gamma =4/3`$, we obtain the simplest perturbative solution given by Equation (5) that is described by the power-law of the time variable. $`\gamma =5/3`$ is given from theoretical argumentation by Buchert and Domínguez budo . They argued kinematic considerations using collisionless Boltzmann equation and derived $`\gamma =5/3`$. The exponent is equivalent to the adiabatic process of an ideal gas. Because we cannot decide on a proportional constant $`\kappa `$ in the equation of state from past discussion, we choose several values. In this paper, instead of $`\kappa `$, we write an initial ($`a=10^3`$, i.e. $`z=1000`$) Jeans wavenumber, given by Equation (6).

Here we show how we set up the initial condition in the pressure model. We adjust the initial condition in the pressure model to be the same as that in ZA: the initial peculiar velocity in the pressure model is same as that in ZA. The procedure for setting up initial condition was shown in our previous papers moritate ; tate02 .

### III.1 Cross-correlation coefficient

First we calculate the cross-correlation coefficient of density fields. The Cross-correlation coefficient was used for the comparison of the resulting density fields cms ; mps ; msw ; Buchert94 ; Melott95 ; Karakatsanis97 . The cross-correlation coefficient is defined by

$$S\frac{\delta _1\delta _2}{\sigma _1\sigma _2},$$
(12)
where $`\sigma _i`$ means the density dispersion of model $`i`$,

$$\sigma _i\sqrt{\delta _i^2}.$$
(13)
$`S=1`$ means that the pattern of density field of two models coincide with each other. In linear regime, the density dispersion remains $`\sigma 1`$. Although we develop the structure until it becomes strongly nonlinear regime ($`\sigma _{Nbody}>1`$), we especially analyze it in the quasi-nonlinear regime ($`\sigma _{Nbody}1`$).

Figure 1-6 shows a comparison of N-body density fields with those predicted by various Lagrangian approximations. First, we notice scale-free spectrum cases (Figure 1-3). As in the past analyses, TZA shows better performance than ZA. Our analyses also show similar tendency, i.e. our analyses do not contradict past analyses.

In the pressure model, the performance strongly depends on the polytrope exponent $`\gamma `$ and Jeans wavenumber. In the case of $`\gamma =4/3`$, when we set the initial Jeans wavenumber to be small, even if in linear regime, the approximation deviates from an N-body simulation. Only for the case of $`K_J=64`$ the approximation shows better performance than ZA in a quasi-nonlinear regime. We notice that the result strongly depends on the Jeans wavenumber in the case of $`\gamma =4/3`$: When we slightly change the value of the Jeans wavenumber, the cross-correlation coefficient changes dramatically. In the case of $`\gamma =5/3`$, although the result depends slightly on initial Jeans wavenumber, the pressure model shows a better performance than ZA in a quasi-nonlinear regime. Furthermore, the pressure model also shows a better performance than TZA. However, when we consider scale-free spectrum models, the model does not have typical physical scale: the model has only box size, grid size, and softening length. The trend of the result was unchanged when we changed the box size of the model and the softening parameter. When we change the number of particles, the result changes. From a comparison of Figure 1 and Figure 3, we can notice that the results depend on the ratio of grid size and initial Jeans wavenumber. In our calculation, we found out that it was good to set up the value of $`\kappa `$ so that the initial ($`z=1000`$) Jeans wavenumber $`K_J`$ was $`N^{1/3}/4K_JN^{1/3}`$. For example, in the case of $`N=64^3`$, as we see in Figure 1 and 2, it is good to choose the initial Jeans wavenumber $`16K_J64`$.

Next we consider SCDM and LCDM models (Figure 4\- 6). In these models, the difference between ZA and TZA becomes very small. Because the initial density spectrum in the CDM model dumps power in a small scale, the cutoff in spectrum weakly affects the formation of caustics, as we saw in the case of $`P(k)k^1`$. From Figure 4, we can see that the effect of pressure improves the approximation in the quasi-nonlinear stage. Also, in both the SCDM and LCDM models, the case of $`\gamma =4/3`$ shows deviation from ZA in linear regime. On the other hand, the case of $`\gamma =5/3`$ shows that the cross-correlation coefficient becomes almost the same in linear regime. In the case of $`\gamma =4/3`$, when we choose a small initial Jeans wavenumber (for example, $`K_J=16`$), although we can improve the approximation much more in quasi-nonlinear stage than in large Jeans wavenumber ($`K_J=32,64`$) cases, the approximation changes slightly for the worse in the linear stage. On the other hand, when we choose $`\gamma =5/3`$, although the effect seems small, we can obtain a improved solution both in linear and in quasi-nonlinear stages.

When the model evolves to strongly nonlinear regime, the trend of the solutions change. In the linear stage, the case of $`\gamma =4/3`$ shows deviation from ZA. However in strongly nonlinear regime, though the Lagrangian approximation generally become worse, the case of $`\gamma =4/3`$ shows a rather good result (Figure 4 (a) and (c)). This tendency was unchanged even though the course-grained length was changed (Figure 5 and 6).

In both the SCDM and LCDM models, when we choose a small initial Jeans wavenumber $`K_J`$, although the approximation is improved after the quasi-nonlinear stage, the reasonable range of Jeans wavenumbers seems wide. The strict limitation to the value of $`\kappa `$ or the initial Jeans wavenumber will be given by other physical considerations or by the high-resolution N-body simulation.

From these results, we find that it is reasonable to choose the polytrope exponent $`\gamma =5/3`$ until quasi-nonlinear regime is reached. These results support the suggestion by Buchert and Domínguez budo  and by Domínguez domi0103 . If we have interest in strongly nonlinear regime, although the Lagrangian approximation generally becomes worse, we had better analyze the case of $`\gamma =4/3`$. In this regime, it is necessary to consider whether that approximation will still be valid. In any case, from the cross-correlation coefficient, we can give limitation in the polytrope exponent $`\gamma `$.

Unfortunately, in these results, we cannot give a strict limit to the proportion coefficient $`\kappa `$ of the equation of state. When we choose $`\gamma =4/3`$, we show that the result strongly depends on $`\kappa `$, and we can notice a strict limitation. However when we choose $`\gamma =5/3`$, we can hardly judge the best value for $`\kappa `$. In our calculation, we found out that it was good to set up the initial ($`a=10^3`$, i.e. $`z=1000`$) Jeans wavenumber $`K_J`$ as $`N^{1/3}/4K_JN^{1/3}`$. From the range of $`K_J`$, we can obtain a reasonable value for $`\kappa `$. If we choose large value for $`\kappa `$, it becomes hard to form nonlinear structure. On the other hand, if we choose small value for $`\kappa `$, the structure becomes almost the same as the structure which was obtained by ZA.

Although the cross-correlation coefficient is one thing which is good for checking the accuracy of the approximation, it is not enough for checking. Now we consider two samples A and B. The density contrast of the samples is given by $`\delta _A`$ and $`\delta _B`$, respectively. We assume that the following proportion relation exists between $`\delta _A`$ and $`\delta _B`$:

$$\delta _A\delta _B.$$
(14)
Even if the density contrast is greatly different, as $`\delta _A`$ and $`\delta _B`$, and one shows a highly nonlinear structure and another remains in the linear regime, the cross-correlation coefficient between $`\delta _A`$ and $`\delta _B`$ becomes $`1`$.

Therefore, we must check the accuracy of the approximation with another property. In the next subsection, we analyze the probability distribution function of density fluctuation.

### III.2 Probability Distribution Function of Density Fluctuation

Here, we compare the probability distribution function (PDF) of density fluctuation. In the Eulerian linear approximation, if initial data is given by a random Gaussian distribution, the PDF of density fluctuations will retain its Gaussianity during evolution. On the other hand, in the Lagrangian approximation, there appears a nonlinear effect. In fact, Kofman et al. Kofman94  shows that the PDF of density fluctuation approaches a log-normal function rather than a Gaussian function in the cases of the Lagrangian approximation and N-body simulation. Padmanabhan and Subramanian PadSub93  also discussed the PDF of density fluctuation with the ZA and found a non-Gaussian distribution.

How will the PDF of density fluctuation change if we take the effect of the velocity dispersion into consideration? Figures 7, 8, and 9 shows the PDF of density fluctuation. As in past works, the PDF of density fluctuation becomes log-normal in form in the N-body simulation. In Figures 7, 8, and 9, the cases of $`\gamma =4/3`$ obviously show the different tendency: in these cases, the effect of pressure suppress the growth of positive fluctuation (Figure 7(b), 8 (c),(d), and 9 (c),(d)). When we also consider the PDF of density fluctuation, we can notice that it is not so good to choose $`\gamma =4/3`$ to examine the growth of structure, though the cross-correlation coefficients well show the trend. On the other hand, the case of $`\gamma =5/3`$ well show the trend in the PDF of density fluctuation. Although the difference of distribution between ZA and the pressure model is still small in quasi-nonlinear regime, the effect of the pressure can promote the evolution of nonlinear structure. Therefore the probability of low and high density regions increases in the case of $`\gamma =5/3`$. Furthermore, according to Figure 8(c), the PDF of density fluctuation in the cases of $`\gamma =5/3`$ show that it is much better than that in the TZA case. Of course when we reach a strongly nonlinear regime, it is necessary to consider whether that approximation is still valid or not.

From both the cross-correlation coefficient and PDF of density fluctuation, we can decide that it is reasonable to choose $`\gamma =5/3`$ as the polytrope exponent of the equation of state. However it is hard to decide the proportional parameter $`\kappa `$. From the results in this paper, we cannot give tight limit to $`\kappa `$. To decide the value of $`\kappa `$, we will analyze a high-resolution N-body simulation or consider other physical processes. For example, we will consider the effect of anisotropic velocity dispersion mtm  or the higher order velocity cumulant.

## IV Discussion and Concluding Remarks

We compared two statistical quantities between an N-body simulation and Lagrangian approximations. In our earlier work moritate ; tate02 , we solved the first-order perturbation equations in the homogeneous and isotropic background, and the second- order ones explicitly for the case $`\gamma =4/3,5/3`$ in Einstein-de Sitter Universe. We showed that the difference between the Lagrangian first-order and second-order approximations becomes small in the case of $`\gamma 4/3`$. Therefore, in this paper we consider only the first-order perturbative solution for the case $`\gamma =4/3,5/3`$. Then we carried out similar calculation with ZA and TZA to examine their difference from the past models.

First, we compared these models by the cross-correlation coefficient of the density field between the N-body simulation and Lagrangian approximations. In scale-free spectrum cases, as well as in the previous analyses, TZA shows a better performance than ZA. In the pressure model, the performance strongly depends on polytrope exponent $`\gamma `$ and Jeans wavenumber. In the case of $`\gamma =4/3`$, when we set that initial Jeans wavenumber to be small, even in linear regime the approximation deviates from the N-body simulation. In the case of $`\gamma =5/3`$, although the result slightly depends on the initial Jeans wavenumber, the pressure model shows a better performance than ZA in quasi-nonlinear regime. Furthermore, the pressure model also shows better performance than TZA. In the SCDM and LCDM models, the case of $`\gamma =4/3`$ shows deviation from ZA in linear regime. On the other hand, the case of $`\gamma =5/3`$ shows that the cross-correlation coefficient becomes almost the same in linear regime. When the model reaches a strongly nonlinear stage, although the Lagrangian approximation generally become worse, the case of $`\gamma =4/3`$ shows a rather good result. Of course, in this regime, it is necessary to consider whether that approximation is still valid.

Second, we analyzed the PDF of density fluctuation. The cases of $`\gamma =4/3`$ obviously show a different tendency until quasi-nonlinear regime is reached: in this case, the effect of pressure suppresses the growth of structure. When we also consider the probability distribution of density, we can see that it is not so good to choose $`\gamma =4/3`$ to examine the growth of structure, although the cross-correlation coefficients perform well. On the other hand, the case $`\gamma =5/3`$ shows good tendency in the PDF of density fluctuation. Although the difference of the PDF of density fluctuation between ZA and the pressure model is still small in quasi-nonlinear regime, the effect of the pressure can promote the evolution of nonlinear structure. The difference between the models of Lagrangian approximation becomes small when we calculate the evolution until strongly nonlinear regime is reached. From analyses of the cross-correlation coefficient of density field and the PDF of density fluctuation, we can decide that it is reasonable to choose $`\gamma =5/3`$ as the polytrope exponent of the equation of state.

In this paper, we changed some values of the Jeans wavenumber $`K_J`$ and undertook analysis. Will there be any relations between the ‘nonlinear wavenumber’ $`k_{NL}`$ in TZA and $`K_J`$ ? The correspondence is as follows. For simplification, we consider the correspondence in case of a scale-free spectrum $`P(k)k^n`$. According to the definition of a ‘nonlinear wavenumber’ in TZA, $`k_{NL}`$ is given from Equation (10). In case of the scale-free spectrum $`P(k)=Ak^n`$, the definition becomes

$$\frac{1}{n+1}a(t)^2Ak_{NL}^{n+1}=1.$$
(15)
From this definition, $`k_{NL}`$ is written as

$$k_{NL}a^{2/(1+n)}.$$
(16)
For example, when we choose $`n=1`$, $`k_{NL}`$ becomes

$$k_{NL}\frac{1}{a}.$$
(17)
On the other hand, the Jeans wavenumber $`K_J`$ in the pressure model is given from Equation (6). When we choose $`\gamma =2`$, $`K_J`$ becomes

$$K_J\frac{1}{a}.$$
(18)
There are some different points to consider when we think about time evolution, although the relation seems to be as described above: First, in TZA, $`k_{NL}`$ affects only the initial spectrum. On the other hand, $`K_J`$ affects the evolution of fluctuation. Second, although $`k_{NL}`$ obviously depends on the initial spectrum, we did not clarify the dependence on the initial condition of $`K_J`$. We think that a consideration of the physical process, which was not considered here, or the analysis of the N-body simulation, is necessary for a decision of $`K_J`$, i.e. $`\kappa `$. We will have to think about the correspondence between the adhesion approximation and the pressure model. Buchert et al. bdp  showed how the viscosity term in the adhesion approximation is generated by a pressure-like force. Domínguez domi00 ; domi0106  discussed spatial coarse graining in a gravitating system and derived an evolution equation of ‘adhesion approximation’. In the pressure model, we showed that the density distribution of the pressure model was similar to that of TZA in the previous paper tate02 . The acute characteristic skeleton structure which appeared in the adhesion approximation could not be seen from the calculations in our previous paper. We will consider the relation between the viscosity term in coarse-grained equations and the pressure term in our model. Then we will analyze the correspondence between the viscosity term in the adhesion approximation and the proportional constant $`\kappa `$ in equation of state in pressure model.

In this paper, we analyzed only density distribution. How will peculiar velocity distribution change with the effect of ’pressure’? In ZA, the peculiar velocity is in proportion to the Lagrangian displacement. Then the growth rate of perturbation is independent of scale. Therefore, although the structure becomes nonlinear regime, if the initial condition is given as Gaussian, the peculiar velocity distribution remains Gaussian all the time Kofman94 . However, in the pressure model, the growth rate of the perturbation depends on the scale. Therefore the peculiar velocity distribution will deviate from Gaussian during evolution. Of course, the peculiar velocity distribution in an N-body simulation becomes non-Gaussian ZQSW94 . Does the effect of the pressure cause the occurrence of the non-Gaussian distribution? We think that the time evolution of the peculiar velocity distribution is one of the more interesting problems.

In our model, we introduce the strong simplification that the velocity dispersion is approximately isotropic, i.e. the stress tensor is diagonal and a pressure-like term bdp . However, in general, the velocity dispersion does not remain isotropic in nonlinear regime. Until when is the assumption to ignore anisotropic velocity dispersion reasonable? Maartens et al. mtm  discussed a relativistic kinetic theory generalization which also incorporates anisotropic velocity dispersion. Then they added these effects to the linear development of density inhomogeneity and found exact solutions for their evolution. In a Newtonian description, though the equations is not generally closed, we will consider anisotropic velocity dispersion and the higher-order velocity cumulant and estimate their effects on the evolution of density inhomogeneity.

###### Acknowledgements.

We would like to thank Thomas Buchert, Kei-ichi Maeda, Masaaki Morita, Momoko Suda, and Hideki Yahagi for useful discussion and comment in the work. For usage of COSMICS and $`P^3M`$ codes, we would like to thank Edmund Bertschinger and Alexander Shirokov. We also thank the referee who gave proper comments and advice for modifications. This work was supported in part by a Waseda University Grant for Special Research Projects (Individual Research 2002A-868 and 2003A-089).


format: text

text:
# Arakelian 564: An XMM-Newton view

## 1 Introduction

Narrow-Line Seyfert 1 galaxies (NLS1s) are a subclass of active galactic nuclei (AGNs) defined by their exceptionally narrow optical permitted emission lines (H$`\beta `$ FWHM $`<2000`$ km s-1) and their generally high Fe II/H$`\beta `$ ratios (e.g., Osterbrock & Pogge 1985; Goodrich 1989), in comparison with “normal” broad-line Seyfert 1 galaxies (BLS1s). NLS1s often show extreme AGN properties; their ultraviolet (UV)/optical emission lines place them at one extreme of the Boroson & Green (1992) primary eigenvector. In the X-ray regime, NLS1s are characterized by rapid and large-amplitude variability (e.g., Boller, Brandt & Fink 1996; Boller et al. 1997); the normalised excess variance of variability is systematically larger in NLS1s than in BLS1s (e.g., Leighly 1999a; Turner et al. 1999; Cancelliere & Comastri 2002), despite both classes having similar X-ray luminosity distributions. The X-ray spectra of NLS1s below $``$ 1–1.5 keV often show a prominent soft excess, modeled either by a steep power law or a thermal component (e.g., Brandt 1999; Comastri 2000). Moreover, their hard ($`>2`$ keV) X-ray spectra are generally characterized by steeper power-law photon indices (e.g., Brandt, Mathur & Elvis 1997; Comastri 2000) than those of BLS1s. A likely explanation for the different properties of NLS1s is that they have relatively low masses for their central black holes and high accretion rates (e.g., Czerny et al. 2001; Boroson 2002; Wang & Netzer 2003). Smaller black hole masses can naturally explain both the narrowness of the optical emission lines, which are generated in gas that has smaller Keplerian velocities, and the extreme X-ray variability, since the primary emission would originate in a smaller region around the central engine. Soft photons from the accretion disc may Compton cool electrons in the corona and cause the steep observed photon indices (Pounds, Done & Osborne 1995; Haardt, Maraschi & Ghisellini 1997). In the case of high accretion rates, the surface of the disc is expected to be ionised (e.g., Matt, Fabian & Ross 1993; Ballantyne, Iwasawa & Fabian 2001). The disc will thus produce ionised Fe K$`\alpha `$ features, as have apparently been observed from some NLS1s (e.g., Comastri et al. 1998; 2001; Vaughan et al. 1999b; Leighly et al. 1999b; Turner et al. 2001a).

Arakelian 564 (hereafter Ark 564; $`z=0.0247`$) is one of the X-ray brightest NLS1s (e.g., Brandt et al. 1994; Vaughan et al. 1999a,b). In 2000 June it was the subject of an intensive multiwavelength monitoring campaign that included simultaneous observations with ASCA (Turner et al. 2001b, hereafter T01; Edelson et al. 2002), RXTE (Pounds et al. 2001), Chandra (Matsumoto, Leighly & Marshall 2001; Marshall 2002), XMM-Newton (this paper), FUSE (Romano et al. 2002), and HST (Collier et al. 2001; Crenshaw et al. 2002). Ark 564 represents a good target for XMM-Newton: the low-energy coverage allows accurate modeling of the soft excess, while the relatively large effective area at high energies allows studies of the previously revealed Fe K$`\alpha `$ emission line (e.g., Comastri et al. 2001). Furthermore, the large XMM-Newton count rate allows the best possible studies of rapid variability.

## 2 Data reduction

Ark 564 was observed by XMM-Newton (Jansen et al. 2001) on 2000 June 17 (rev. 96) and 2001 June 9 (rev. 275).111Ark 564 was observed twice because of operational difficulties that prevented the entire accepted exposure from being obtained in one observation.

At the beginnings of both observations the European Photon Imaging Camera (EPIC) instruments \[pn (Strüder et al. 2001) and MOS (Turner et al. 2001)\] were closed, while the Reflection Grating Spectrometer (RGS; den Herder et al. 2001) and the Optical Monitor (OM; Mason et al. 2001) were operated. The EPIC cameras were operated in “Small-Window” mode to mitigate “pile-up” of photons from the source; EPIC data were acquired using the medium optical blocking filter. The data were reduced with the latest version of the XMM-Newton Science Analysis Software (sas; version 5.4.1) using the latest calibration products. The data were filtered to avoid background flares (due to $``$ 100 keV protons; e.g., De Luca & Molendi 2002) which affect large intervals ($``$ 50%) of our observations. For the EPIC data, the tasks epproc and emproc were used to generate valid photon list. Both single-pixel and double-pixel events (patterns 0–4) were used when extracting the pn counts, while patterns 0–12 were used for the MOS. Source counts were extracted in the 0.3–12 keV band from circular regions of radius 90′′ (pn) and 50′′ (MOS); background counts were extracted from off-source regions of radius 45′′ (pn) and 150′′ (MOS). The statistics are clearly dominated by the source counts in any chosen EPIC sub-band. Response functions for spectral fitting to the EPIC data were generated using the tasks rmfgen and arfgen. First-order RGS spectra were extracted using the rgsproc task, which also produced the appropriate response matrices. A summary of the XMM-Newton observations of Ark 564 along with the source counts are presented in Table 1. Note the large numbers of counts obtained with the present XMM-Newton observations of Ark 564; for comparison, each ASCA Solid-state Imaging Spectrometer gathered $``$ 2$`\times 10^6`$ source counts over the whole 35-day observation (T01) and $``$ 12000 in the same time interval as the XMM-Newton observation (see $`\mathrm{\S }`$3.1).

## 3 Spectral analysis

The source counts were grouped into spectra such that each spectral bin contained at least 20 counts to allow $`\chi ^2`$ fitting. X-ray spectra were fitted using the xspec package (version 11.2; Arnaud 1996). The quoted errors on derived model parameters correspond to the 90 per cent confidence level for one interesting parameter (i.e., $`\mathrm{\Delta }\chi ^2=2.71`$; Avni 1976) unless otherwise stated. All spectral fits include absorption due to the line-of-sight Galactic column density of $`N_\mathrm{H}=6.4\times 10^{20}`$ cm-2 (Dickey & Lockman 1990). Hereafter we adopt $`H_0`$=70 km s-1 Mpc-1 in a $`\mathrm{\Lambda }`$-cosmology with $`\mathrm{\Omega }_\mathrm{M}`$=0.3 and $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7.

At first the MOS1 and MOS2 spectra were fitted separately to check for cross-calibration uncertainties; these have been found to be $``$ 3% (Kirsch 2003). We found generally good agreement between MOS1 and MOS2, comparable with the value quoted above. Therefore in the following we will present the spectral results of both cameras referred to as MOS.

Cross-calibration uncertainties between pn and MOS are generally at most $``$ 13% (Kirsch 2003). However, the pn and MOS spectra of Ark 564 in both observations clearly differ at energies below $``$ 0.6 keV; this is probably due to residual calibration problems in the “Small-Window” mode. A similar difference in the pn/MOS behaviour at low energies has been found for PG 1211$`+`$143 (see Fig. 3 of Pounds et al. 2003). Therefore, in the following the pn and MOS data will be fitted together with the same model in the 0.6–10 keV band, leaving the relative normalisations free to vary. We note, however, that some discrepancies between the pn and MOS seem to be present up to $``$ 0.9 keV. The extension of our analysis down to 0.3 keV has the advantage of providing a slightly better constraint on the temperature of the soft thermal component without changing our main results.

### 3.1 Spectral results: The 2000 June observation

A single power-law model provides a poor fit to the 0.6–10 keV pn$`+`$MOS data; strong residuals at low energies are evident in both instruments, as shown in Fig. 1.

The soft excess of Ark 564 rises sharply below $``$ 1.5 keV similarly to the BeppoSAX observation (see Fig. 6 of Comastri et al. 2001). The presence of a soft component, parameterised either by a blackbody or a power law, is common across the NLS1 population (e.g., Comastri 2000; Vaughan et al. 2002). If thermal, this component is usually ascribed to the high-energy tail of the accretion disc responsible for the UV bump detected in the majority of NLS1s (see, e.g., Soria & Puchnarewicz 2002 for discussion), although the rapid variability of the soft component detected from some NLS1s casts doubts on simple versions of the accretion-disc interpretation. Fitting the soft excess with a power law is inadequate; there is clear continuum spectral curvature that is not accounted for by a power law. Spectral curvature with similar shape has been detected in some other Seyfert galaxies (see, e.g., Pounds & Reeves 2002 for a review; Collinge et al. 2001) and also in the narrow emission-line quasar PG 1211$`+`$143 (Pounds et al. 2003). The addition of a blackbody component is preferred by the EPIC data and finds support from previous BeppoSAX results (Comastri et al. 2001). This component has $`kT=138_3^{+2}`$ eV \[see model (1) in Table 2\].

At higher energies, a power law with photon index $`\mathrm{\Gamma }=2.55\pm 0.02`$ dominates the spectrum. This model leaves significant residuals at soft X-ray energies. The shape of these residuals (see Fig. 2) indicates a count deficit at $``$ 0.7 keV.

We therefore added an absorption edge to the previous model \[see model (2) in Table 2\]. The edge provides a highly significant improvement in the fit quality ($`\mathrm{\Delta }\chi ^2`$/$`\mathrm{\Delta }`$d.o.f.=163/2). Its rest-frame threshold energy, 731$`\pm 11`$ eV, is consistent with the edge being due to vii (which has an energy of 739 eV). The 68, 90, and 99 per cent confidence contours for the edge energy vs. optical depth (in the observed frame) are shown in Fig. 3. The edge is still required when other continuum models for the soft component (e.g., a power law, a double blackbody) are chosen, although the fit is significantly worse.

The edge is also present in the Chandra High-Energy Transmission Grating (HETG) data taken simultaneously with the XMM-Newton observation (see the spectrum shown in Fig. 3 of Matsumoto et al. 2001). The edge rest-frame threshold energy in the Chandra observation (712$`{}_{4}{}^{}{}_{}{}^{+3}`$ eV; Matsumoto et al. 2001; Marshall 2002) is not in agreement with the nominal energy of the vii K-edge. The hypothesis that Chandra observed a highly redshifted vii K-edge with a velocity of $``$ 1000 km s-1 appears unlikely since the vii absorption lines detected by Chandra have no measurable velocity shifts (Matsumoto et al. 2001). We have performed an analysis of the RGS data searching for the spectral features reported by Matsumoto et al. (2001) based upon HETG data. However, we do not find any clear absorption lines (at most $`\mathrm{\Delta }\chi ^2`$/$`\mathrm{\Delta }`$d.o.f.=10/2 for a narrow absorption line) despite the high counting statistics of the RGS data. This is presumably due to the lower spectral resolution of the RGS; the HETG spectral features are unresolved even with the higher spectral resolution of the HETG. Therefore, below we use the RGS mainly to check the results derived from the EPIC spectral fitting rather than presenting a detailed analysis of the RGS data.

The edge is also present in the RGS spectra, although at a much lower significance level; spectral simulations suggest that this result is due to the combination of lower statistics and the presence of chip gaps in the RGS. The underlying continuum measured by the HETG is well parameterised by a model similar to the one used for our XMM-Newton observation, even though the HETG bandpass is significantly narrower than the XMM-Newton one. While the HETG best-fitting photon index ($`\mathrm{\Gamma }=2.56\pm 0.06`$) is consistent with the one measured by XMM-Newton, the blackbody component has a lower temperature ($`kT=124\pm 3`$ eV). The different blackbody temperature may be the cause of the different edge energy in the HETG observation.

Adopting the blackbody plus power-law model and adding a warm absorber (absori in xspec; Magdziarz & Zdziarski 1995) to fit the absorption feature of Ark 564 does not provide good results; significant residuals are still present at the energy corresponding to the vii edge. The ionisation parameter, defined as $`\xi =\frac{L}{nR^2}`$, where $`L`$ is the source luminosity, $`n`$ the warm absorber density, and $`R`$ the distance of the absorber from the photoionising source, is 73$`{}_{35}{}^{}{}_{}{}^{+43}`$ erg cm s-1. For reasonable assumptions of the physical parameters, the strongest bound-free soft X-ray absorption features in the warm absorber model are the vii and viii K-edges (e.g., Reynolds & Fabian 1995; Reynolds 1997; George et al. 1998). The relatively large optical depth of the detected vii edge probably causes the failure of the warm absorber model, since absori also predicts the presence of the viii edge which is not observed in the data ($`\tau `$$`<`$0.08). Although data with better resolution and higher signal-to-noise ratio are required to provide a physical justification for the presence of the vii edge and the absence or weakness of the viii edge, we note that similar spectral results have been found for the Seyfert 1 galaxy H 1419$`+`$480 using XMM-Newton data (Barcons, Carrera & Ceballos 2003). Although the different spectral resolution does not allow a proper comparison of the absorption features at different wavelengths, there are claims that the UV and X-ray absorbers in Ark 564 are physically related and possibly identical (e.g., Romano et al. 2002; Crenshaw et al. 1999, 2002; Matsumoto, Leighly & Marshall 2002). However, the $`N_{\text{}\text{vii}}`$ value we derive from the observed K-edge optical depth ($`1.4\times 10^{18}`$ cm-2) is about an order of magnitude larger than the value predicted in Table 4 of Crenshaw et al. (2002) using photoionisation models and assuming the UV absorber to be a single zone. It is likely that a more complex structure for the absorber is required to reproduce both the UV and X-ray data.

In 2000 June Ark 564 was also observed by ASCA (T01) with a 35-day observation (from 2000 June 1 to July 6). Unfortunately, the lower statistics in the portion of ASCA data (re-analysed using standard procedures; $``$ 12000 counts per instrument) taken simultaneously with our XMM-Newton observation and the exclusion of the data below 0.7 keV do not allow one to reveal the vii edge detected by XMM-Newton.222See $`\mathrm{\S }`$2 of T01 for a discussion of the ASCA calibration problems that led the authors to avoid the spectral band below 0.7 keV. However, the spectral parameters for the blackbody ($`kT=146\pm 9`$ eV) and the power-law ($`\mathrm{\Gamma }=2.46\pm 0.03`$) components are in relatively good agreement with our results.

At high X-ray energies, the long ASCA observation indicates that the iron K$`\alpha `$ emission line is complex, best fitted using a relativistic iron line produced by an accretion disc in either the Schwarzschild or Kerr metrics (xspec models diskline and laor, respectively; Fabian et al. 1989; Laor 1991), and variable on time scales of a few days. The analysis of the ASCA data exactly simultaneous with our XMM-Newton observation, however, indicates only marginal evidence ($`\mathrm{\Delta }\chi ^2`$/$`\mathrm{\Delta }`$d.o.f.=6/2) for an ionised Fe K$`\alpha `$ line (EW=113$`{}_{79}{}^{}{}_{}{}^{+181}`$ eV). This is possibly caused by the relatively poor statistics in the short time interval simultaneous with the XMM-Newton observation.

The addition of a narrow ($`\sigma =10`$ eV) emission line improves the EPIC fit quality by $`\mathrm{\Delta }\chi ^2`$=8 at most when all of the line parameters, except for its normalisation, are frozen. We obtain EW=32$`{}_{17}{}^{}{}_{}{}^{+46}`$ eV, EW=35$`{}_{31}{}^{}{}_{}{}^{+30}`$ eV, and EW=43$`{}_{35}{}^{}{}_{}{}^{+36}`$ eV for the neutral (6.40 keV), He-like (6.70 keV), and H-like (6.97 keV) Fe K$`\alpha `$ lines, respectively. We tried to account for the small residuals at high energies, having a “bump-like” shape (see Fig. 4), using a relativistic iron line produced by an accretion disc. The diskline model improves the fit by $`\mathrm{\Delta }\chi ^2`$/$`\mathrm{\Delta }`$d.o.f.=12/2 when all the line parameters are frozen to the default values333We assumed the line originates between 10 and 1000 gravitational radii and adopted an emissivity law $`r^q`$ with $`q=2`$ for the illumination pattern of the accretion disc; the inclination angle of the disc was fixed to 30$`\mathrm{°}`$. excluding the line energy and its normalisation, which are left free to vary. The line rest-frame energy is 6.10$`{}_{0.17}{}^{}{}_{}{}^{+0.31}`$ keV, while its intensity is EW=98$`{}_{46}{}^{}{}_{}{}^{+48}`$ eV. Both the line energy and EW are lower than the values found by T01. This is probably due to the fact that this model fails in reproducing the peak of the residuals at $``$ 7 keV in the present XMM-Newton data. Using the laor model444We fixed the innermost radius as the last stable orbit for a Kerr black hole and the outer radius at the maximum value allowed by the model (400 gravitational radii). We fixed the emissivity index to 3 and the inclination angle of the disk to 30$`\mathrm{°}`$. the improvement is slightly more significant ($`\mathrm{\Delta }\chi ^2`$/$`\mathrm{\Delta }`$d.o.f.=18/2). The line rest-frame energy, considerably higher (7.07$`{}_{0.20}{}^{}{}_{}{}^{+0.27}`$ keV) than in the previous model, and its EW (471$`{}_{86}{}^{}{}_{}{}^{+175}`$ eV) are consistent with the best-fit model of T01, where the line peaks at $``$ 7 keV. However, in contrast with the ASCA results shown by T01, the presence of a complex Fe K$`\alpha `$ line is not strongly motivated by the present XMM-Newton data, either by the residuals to a power-law fit or by the improvement in the $`\chi ^2`$ after the inclusion of this feature. It is possible that residual uncertainties either in the calibration of XMM-Newton spectra at high energies in “Small-Window” mode or in ASCA, coupled with a rapid change in the intensity and shape of the line (already probed by the longer ASCA observation), are responsible for the different results we obtained. Further investigations on this issue are required.

When a reflection component is added to the blackbody plus power-law model, a much steeper photon index is obtained (as expected, given the known correlation of these parameters in the reflection models available in xspec; see Vaughan & Edelson 2001 and Perola et al. 2002 for discussion), but the fit does not improve significantly. The generally flatter photon index in MOS1 data (see Molendi & Sembay 2003) and, possibly, problems in the pn-MOS cross-calibration using the “Small-Window” mode could produce a spurious reflection component in the MOS data. The significantly higher statistics collected by the EPIC pn at high energies reveal no evident additional spectral components and indicate that model (2) in Table 2 is a good parameterisation of the Ark 564 spectrum in the 2000 June observation. However, motivated by the discoveries of sharp drops in the XMM-Newton spectra of two other NLS1s, 1H 0707$``$495 (Boller et al. 2002) and IRAS 13224$``$3809 (Boller et al. 2003), we fitted the 2–10 keV data using the best-fitting power law obtained in the 2–7 keV band (thus avoiding any problems with the vii edge and the soft thermal component). In contrast with the results obtained for the other two NLS1s, we find only marginal evidence for a drop at $``$ 7–8 keV ($`\mathrm{\Delta }\chi ^2`$/$`\mathrm{\Delta }`$d.o.f.=6/2 when the drop is parameterised by an absorption edge). This is a further indication that the present XMM-Newton data do not require more complex modeling at high energies.

### 3.2 Spectral results: The 2001 June observation

The 2001 June XMM-Newton EPIC observation of Ark 564 is characterized by lower statistics (see Table 1) than the 2000 June observation. The exposure time, after the screening procedure, is relatively low, and the source has a lower flux level (see $`\mathrm{\S }`$4). Following $`\mathrm{\S }`$3.1, a blackbody plus power law is used to model the broad-band X-ray continuum of Ark 564; the blackbody temperature and the power-law photon index \[see model (1) in Table 3\] are consistent with the corresponding model in the 2000 June observation (see Table 2). The poor quality of the fit, however, is suggestive of additional spectral complexity, in particular below $``$ 1 keV.

Based on the 2000 observation and the shape of the residuals in the 2001 observation, we tried to account for the poor spectral fit by adding an absorption edge \[model (2) in Table 3\].

The edge is significant statistically ($`\mathrm{\Delta }\chi ^2`$/$`\mathrm{\Delta }`$d.o.f.=68/2); its rest-frame threshold energy, $`727_{16}^{+15}`$ eV, is consistent with the vii K-edge (see Fig. 5), as in the 2000 June observation, and also the edge optical depth agrees with the previous EPIC observation. Some indications of this spectral feature are also present in the RGS data. There is no clear evidence for additional spectral features in the 2001 XMM-Newton data. An upper limit of $`\tau =0.12`$ is derived for the viii optical depth; adopting a warm absorber model to reproduce the vii edge provides a worse spectral fit and a poorly constrained ionisation parameter.

Although in the 2001 June observation of Ark 564 the hard X-ray photon statistics are limited, we can place some constraints on the presence of Fe K$`\alpha `$ emission lines (both neutral and ionised). We find marginal evidence ($`\mathrm{\Delta }\chi ^2`$/$`\mathrm{\Delta }`$d.o.f.=12/2; the width of the line is frozen to $`\sigma =10`$ eV) for a Gaussian Fe K$`\alpha `$ line. The line rest-frame energy ($`6.16_{0.02}^{+0.05}`$ keV), lower than that expected for the neutral Fe K$`\alpha `$ emission line, coupled with its low significance casts some doubts on the nature of this feature. If we freeze the line energy at 6.40 keV, 6.70 keV, and 6.97 keV for the neutral, He-like, and H-like Fe K$`\alpha `$ lines, we obtain EW$`<71`$ eV, EW=70$`{}_{61}{}^{}{}_{}{}^{+60}`$ eV, and EW$`<114`$ eV, respectively, but none of these features is significant statistically. The absence or weakness of iron lines is supported by the lack of a strong reflection component, as the relatively good fit obtained at high energies with model (2) also suggests.

## 4 Flux variability

Between the 2000 and 2001 observations, Ark 564 shows significant flux variability. In the 2000 June observation, the observed soft ($`F_{0.52\mathrm{keV}}`$) and hard ($`F_{210\mathrm{keV}}`$) band fluxes are (4.7–5.4)$`\times 10^{11}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ and (2.1–2.4)$`\times 10^{11}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, respectively, depending on the EPIC camera used and adopting the best-fitting models reported in Table 2. These correspond to Galactic absorption-corrected luminosities of (8.8–10.0)$`\times 10^{43}`$ erg s-1 and (3.0–3.4)$`\times 10^{43}`$ erg s-1, respectively. In 2001 June the source shows a flux level lower by $``$ 40–50 per cent. We found $`F_{0.52\mathrm{keV}}`$=(3.0–3.2)$`\times 10^{11}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ and $`F_{210\mathrm{keV}}`$=(1.1–1.2)$`\times 10^{11}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ \[model (2) in Table 3\], corresponding to intrinsic luminosities of (5.7–6.0)$`\times 10^{43}`$ erg s-1 in the soft band and (1.7–1.8)$`\times 10^{43}`$ erg s-1 in the hard band. In both observations, the presence of a larger luminosity in the soft band than in the hard band argues against a simple model where the soft excess is just due to reprocessed hard X-ray emission.

For each observation, pn light curves were extracted in the 0.2–10 keV range for continuous periods of 9.9 and 11.5 ks from the 2000 and 2001 data sets, respectively, and binned to 50-s time resolution (see Fig. 6, top panels).555We used the total (whole interval) light curves irrespective of the periods of “flaring background” because the source count rate is much higher than the background even during the flares. The background flares are weak compared to many seen by XMM-Newton, and the exclusion of these intervals for the spectral analysis was conservative. Their effect on the PSD derived in $`\mathrm{\S }`$5 is negligible. Similarly, the exclusion of the 0.2–0.3 keV data from the spectral analysis was a conservative choice. The light curves were background subtracted and exposure corrected (i.e., corrected for telemetry drop outs). Analyses of the pn light curves in different energy bands (0.2–10 keV, 0.2–0.7 keV, 0.7–2 keV, and 2–10 keV; see Fig. 6) show that the source is highly variable, both at low and high energies. The rapid variability of the soft component and its relatively high temperature (see $`\mathrm{\S }`$3.1 and 3.2) suggest that it is unlikely to be simple thermal radiation from an accretion disc. It is possible that it is quasi-thermal emission originally from the accretion disc that has been upscattered in a hot surface layer of the disc or the corona. Alternatively, the soft component might be due to turbulent Comptonisation occurring close to the black hole (Socrates, Davis & Blaes 2003).

## 5 Power spectral density function

The high-frequency power spectral density function (PSD) of Ark 564 was estimated from the pn light curves using the method of Vaughan, Fabian & Nandra (2003). pn 2000 and 2001 periodograms were computed using standard Fourier methods. The two periodograms appeared consistent with one another, and therefore a single binned periodogram was produced by combining the two periodograms and binning, using the method of Papadakis & Lawrence (1993), such that there are 15 periodogram estimates per frequency bin. Figure 7 shows the average periodogram. At high frequencies the power is dominated by the Poisson noise.

At low frequencies the source variations dominate with a steep, approximately power-law spectrum rising toward low frequencies.

From previous monitoring with RXTE and ASCA it is known that the low-frequency PSD is well represented by a broken power law, with a power-law index of $`0`$ below a frequency $`f_{\mathrm{br}}1.6\times 10^6`$ Hz and an index of $`1.2`$ above the break (Markowitz et al. 2003). The XMM-Newton data at much higher frequencies show a much steeper spectrum and imply the PSD must show a high-frequency break as well, confirming the suggestion of Papadakis et al. (2002) on the basis of ASCA data. In order to estimate the high-frequency PSD, these data were fitted using the Monte Carlo procedure of Vaughan et al. (2003) with a doubly broken power-law model. The power spectral index was fixed at $`0`$ below $`10^6`$ Hz and $`1.2`$ above this frequency. Then a second, high-frequency break was included in the model with both the break frequency ($`f_{\mathrm{br}}`$) and the spectral index above the break ($`\alpha `$) as free parameters in the fitting.666The low-frequency PSD needed to be specified in the model as the fitting procedure accounts for biases in the periodogram, notably “red-noise leak”, that depend on the PSD at low frequencies. See Vaughan et al. (2003) for details on this.

This model provided a statistically acceptable fit to the data ($`\chi _\nu ^2<1.0`$), but the free parameters were poorly constrained. The PSD is steep, with an index $`\alpha >1.95`$, but the location of the break to this steep slope is only constrained to be $`f_{\mathrm{br}}<8\times 10^4`$ Hz (90 per cent confidence limit for one interesting parameter). This fit thus supports the claim by Papadakis et al. (2002) that Ark 564 shows a doubly broken PSD, although the inferred position of the break is a factor $``$ 2 lower than that obtained in their analysis.

Curiously, the doubly broken PSD resembles that seen in Cyg X-1 in its low/hard state, contrary to the expectation that Ark 564 might look more like the high/soft state (which only shows the high-frequency break). However, it is vital to constrain the high-frequency PSD and to confirm the second break in the PSD before a fair comparison can be made between Ark 564 and Galactic black hole candidates (which can also show “very high” and “intermediate” state PSDs).

## 6 Summary

We have reported two XMM-Newton observations of the bright NLS1 Ark 564 taken one year apart. The principal results are the following:

* A blackbody with $`kT`$ 140–150 eV and a power law with $`\mathrm{\Gamma }`$ 2.50–2.55 provide a good parameterisation of the X-ray continuum. No significant spectral changes are observed between the two observations, although the X-ray flux in the second observation is lower by $``$ 40–50 per cent.
* In both observations a significant absorption edge at $``$ 0.73 keV, corresponding to the vii K-edge, is detected.
* We find tentative evidence for an iron emission line (fitted using a laor model) in the 2000 June observation.
* The present XMM-Newton data support the idea that the PSD shows two breaks. A longer observation is required to place stronger constraints on the location of the high-frequency break.

## Acknowledgments

The work reported herein is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). We gratefully acknowledge the financial support of NASA LTSA grant NAG5-13035 (CV, WNB) and NASA grant NAG5-9939 (CV, WNB). CV also acknowledges partial support from the Italian Space Agency under contract ASI I/R/073/01. ACF thanks the Royal Society for support. We thank D. Alexander, A. Comastri and G. Matt for interesting suggestions, S. Molendi and P. Ranalli for useful discussions about XMM-Newton calibrations, L. Angeretti for help with the SM macros, and the referee for her/his thoughtful comments and suggestions which improved the quality of the paper.


text:
# Efficiency Gains from Time Refinement on AMR Meshes and Explicit Timestepping

## 1 Introduction

### 1.1 Block-Structured AMR

Adaptive mesh refinement on rectangular grids (henceforth AMR) was introduced in bergeroliger , and improved for conservation laws in bergercollela , henceforth BC89. In the patch-based meshes of the sort described in BC89, the patches increase in resolution by a fixed even integer factor $`N`$. One can place a finer patch anywhere in the domain of a ‘parent’ patch of one fewer level of refinement. A patch is not required to have only a single parent, but must be completely contained within patches of the next lowest level of refinement. Note that these meshes are non-conforming; the face of a zone in a parent patch will abut $`N`$ faces in the child patch. A final restriction in the nesting of the meshes is that there must be at least one zone of the next lower level refinement about the perimeter of a patch.

Another mesh we will consider here is an oct-tree mesh (quad-tree in 2-d, binary tree in 1d), such as is implemented in the Paramesh package paramesh  used in the Flash code flashcode . This oct-tree mesh is a more restrictive version of an $`N=2`$ patch-based mesh as described in BC89. If a block needs additional resolution somewhere in its domain, the entire block is halved in each coordinate direction, creating $`2^d`$ children, where $`d`$ is the dimensionality. Leaf blocks are defined to be those blocks with no children, and are thus at the bottom of the tree — they are the finest-resolved blocks in their region of the domain. Frequently, only leaf blocks are evolved to compute the solution to the equations, since a refined parent block’s domain is completely spanned by its children.

The only difference between the two meshing approaches of immediate interest is the resulting different refinement patterns. We will use ‘patch’ and ‘block’ interchangeably in this paper.

### 1.2 Time Refinement

In BC89, the timestep set by the data on the finest mesh is used to evolve that data, and data on the coarser meshes is evolved at a multiple thereof so that there is a constant ratio at each level $`l`$ of $`\mathrm{\Delta }t_l`$ to $`\mathrm{\Delta }x_l`$. The assumption here is that there is one roughly spatially constant characteristic speed throughout the entire domain, so that the maximum allowable timestep at any given resolution is directly proportional to the size of the mesh for any given block or patch. When coupled with the assumption in structured AMR of some fixed jump in refinement between levels, this makes for a very natural time evolution algorithm, shown pictured in Figure 1 for a mesh with three different levels of refinement, with resolution jumps by constant factors of $`N`$; shown is $`N=2`$.

Here the largest blocks are evolved at some system timestep $`dt`$, and smaller blocks are ‘subcycled’ at proportionally smaller timesteps. This defines a ‘work function’ for each block; the finest blocks must be evolved every sub-timestep so we take their work value to be 1 times the number of zones in the block or patch; the blocks one level of refinement ‘up’ need only be evolved every $`N`$ sub-timesteps, so that their work value is $`1/N`$ times the number of zones, etc. The work function for an entire mesh is the sum of the work values of each block or patch in the mesh.

There are costs associated with this time refinement (hereafter TR). Memory is needed to store information at multiple timesteps. There are overheads from extra copies and time-centering of fluxes. The modified time-structure of work leads to load-balance issues in parallel jobs. Further complicating parallel performance is increased communication complexity (although, it is to be pointed out, not necessarily increased communication).

Nonetheless, one might hope that these costs are outweighed by the time savings of not evolving large blocks at unnecessarily small timesteps; in the example of Figure 1, of evolving the larger blocks at timesteps of $`dt`$ or $`dt/2`$ instead of $`dt/4`$. As a first step to quantify the possible benefits, we estimate the reduction in computational cost in simple cases §2. We then use the same approach to examine meshes from simulations performed with a tree-based mesh in §3. In our final section we summarize our results.

## 2 Simple Mesh Configurations

Here we calculate both the number of evolved blocks in a simple mesh, and a weighted sum representing the ideal amount of work done by a TR method, using the work function described in the previous section. We then calculate a work ratio, $`R`$ — the amount of work that would be done by the idealized TR divided by that done with no time refinement. With no time refinement, each block must be stepped through each sub-timestep, so that the amount of work done is simply the number of blocks; thus, the work ratio is simply (TR work function)/(number of blocks). For $`R=1`$, there is no reduction in work; for $`R<1`$, TR reduces the amount of computational work.

One can interpret the work ratios as performance metrics for the TR, assuming that – all physics benifits from the time subcycling in proportion to the reduction in number of blocks evolved each step; the memory overhead from TR is unimportant; all larger blocks actually *can* be evolved at timesteps of larger size in proportion to their physical size; there is no single-processor overhead from TR from memory copies or flux averaging; there is no parallel overhead from increased complexity in communications; and there is no parallel from increased load-balancing issues.

### 2.1 Point refinement

The best case for efficiency gains for spatial refinement is clearly one isolated point of refinement. For a patched-based mesh, we imagine refinements as shown on the left of in Figure 2.

We begin with domains of length one in all directions. The completely unrefined domain is defined to be at level $`l=1`$ of refinement. Consider placing increasingly fine patches at the corner, until we resolved the finest scale $`\mathrm{\Delta }x`$ we wished. If this requires $`L1`$ more levels of refinement, each decreasing the zone size by an integer factor $`N`$, then we have $`\mathrm{\Delta }x(1/N)^{L1}`$. We will assume $`\mathrm{\Delta }x1/2`$.

We consider the mesh in terms of the smallest uniform unit — for the oct-tree mesh, this is a single block, which will be of size $`n_x\times n_y\times n_z`$ zones. For the patch-based mesh, since the patches can be of arbitrary size (and shape), we consider zones individually. (Because we are not modelling guardcell filling, we can safely ignore the fact that these zones are actually components of patches). Thus, in the results given below, an oct-tree mesh with (say) $`8\times 8`$-zone blocks at a maximum refinement $`L=5`$ has the same resolution as a patch-based mesh with $`L=8`$.

The amount of work required by a non-TR code with only explicit or local solves will, by assumption, be the same for each block, so that $`W_{\mathrm{noTR}}=N_{\mathrm{blocks}}`$. The amount of work with time refinement, $`W_{\mathrm{TR}}`$, will be a weighted sum of blocks. For the pointwise-refined patch mesh, the number of blocks will simply be $`N_{\mathrm{blocks}}=L`$, as there is only one block per level. The amount of TR work is

$$W_{\mathrm{TR}}=\underset{l=1}{\overset{L}{}}1\left(\frac{1}{N}\right)^{Ll}\frac{N}{N1}.$$
(1)
Thus the work ratio will be

$$R=\frac{W_{\mathrm{TR}}}{W_{\mathrm{noTR}}}=\frac{W_{\mathrm{TR}}}{N_{\mathrm{blocks}}}=\frac{N}{L(N1)}.$$
(2)
For ideal spatial AMR, where one can do all the refinement with only one jump, $`L=2`$, and so the amount of work done by a TR algorithm is bounded from below at $`1/2`$ of the non-TR work. At the other limit, for a much less aggressive AMR with $`N=2`$, then the work can be made an arbitrarily small fraction of the non-TR algorithm, with $`R=2/L`$ — but note that this work ratio is achieved only by operating on $`L/2`$ times as many blocks as in the best case for spatial AMR.

The oct-tree meshes refining on a point is shown on the right of Figure 2. In this case, there are $`2^d`$ highest refined blocks in the corner, with the rest of the $`2^d1`$ surrounding blocks at the next highest refinement, surrounded by the $`2^d1`$ surrounding blocks at the next highest level of refinement, and so on.

Thus the total number of leaf blocks is

$$N_{\mathrm{blocks}}=(2^d)+\underset{l=L1}{\overset{1}{}}\left(2^d1\right)=2^d(L1)L+2$$
(3)
Weighting them by the amount of work,

$$W_{\mathrm{TR}}=(2^d)+\underset{l=L1}{\overset{1}{}}\left(\frac{(2^d1)}{2^{(Ll)}}\right)2^{(d+1)}1$$
(4)
making the work ratio

$$R=\{\begin{array}{cc}3/L& 1d\hfill \\ 7/(3L2)& 2d\hfill \\ 15/(7L6)& 3d\hfill \end{array}$$
(5)
As with the patch-based result, this ratio goes to zero for arbitrarily large $`L`$. These results are similar to the $`N=2`$ patch-based result, but TR performs better here, and the spatial refinement worse — both of these are due to the fact that the oct-tree mesh generates more intermediate-level blocks.

### 2.2 Planar Interface Refinement

The refinement of an interface is shown on the right of Figure 2. In the patch-based case, we continually place a grid of $`N`$-by-1 (in 2d) or $`N^2`$-by-1 (in 3d) patches along the interface, until the required resolution is achieved.

In this case, performing the same calculation as in the previous section, one obtains

$$R1\frac{N1}{N^d1}.$$
(6)
Here, there is a fixed lower bound for the amount of work the TR can achieve. In the spatially-optimal large-$`N`$ limit, no work is saved at all: $`R1`$. At the other limit, for $`N=2`$, in 2d, $`R2/3`$; in 3d, $`R6/7`$.

In the oct-tree mesh we begin with one block at the coarsest level. It must be divided into 4 in this 2D example, or, in general, $`2^d`$. Half of these blocks will be further refined. This continues until we reach the maximum level of refinement. The work ratio one finds is

$$R=\{\begin{array}{cc}7/9& 2d\hfill \\ 45/49& 3d\hfill \end{array}$$
(7)
In the point-refinement case of the previous subsection, a point of zero volume needed to be refined; as a result, there were the same number of blocks at each level, and thus a significant time savings could be obtained by doing less work at the coarser blocks. However, as we begin to see here, as soon as a non-trivial volume of the mesh needs to be refined, there is significantly less savings to be had.

### 2.3 Circular Region Refinement

The loss of efficiency gains when a non-zero fraction of the mesh must be refined is even clearer when a region, rather than an interface, is fully refined. In Fig. 3 we see the results of fully refining the interior of a quarter-circle with the center at one of the corners of the domain. Clearly, the number of finest blocks greatly outnumber intermediate or large blocks, so one might guess that there is very little efficiency gain that can be had from reducing work on the larger blocks.

Because in this case the refinement pattern is complicated enough that the process must be iterated to check that each zones neighbors are no further than one level of refinement appart, we do not provide analytic work ratios. Tables 1 and 2 show the work ratios for an Oct-Tree mesh and an $`N=2`$ patch-based mesh in refining a circular region of radius $`r`$. Again, the $`r=0`$ results reproduce the expected point refinement, but as soon as a non-zero radius must be refined, the efficiency gains drop significantly further than in the case of only refining an interface, as more small blocks are needed to refine a region than the interface. In Table 2 we also show results for the patch based mesh with $`N=4`$; we see as in previous sections that for the same resolution, increasing $`N`$ (which increases the spatial efficiency of AMR) decreases the possible gains from time subcycling.

## 3 Meshes from simulations

The calculations of the previous section are for very simple refinement geometries. In this section, we apply the same work function used in §2 to the output of previous actual AMR simulations which use oct-tree based meshes for AMR. We continue to assume the same idealized performance results of the previous section.

We begin with examining results from a standard test problem, a Sedov explosion sedov , as included with the Flash code and described in flashcode . In this simulation, a high pressure at a point causes a spherical shock wave to expand outwards; this is analogous to the circular region analysis of the previous section. The adaptive mesh for different stages of this simulation in 2d are shown in 4.

Results from the meshes shown are tabulated in Table 3. The number of blocks listed in the table is the number of ‘leaf’ blocks – *e.g.*, the blocks that are actually evolved. Also given in the table is the work ratio ($`R`$) and the work ratio of spatial AMR to a uniform mesh at the highest resolution ($`R_{\mathrm{AMR}}=W_{\mathrm{noTR}}/W_{\mathrm{uniform}}`$). We include $`R_{\mathrm{AMR}}`$ to compare the relative importance of performance gains for the spatial refinement and the time subcycling.

TR provides a large performance gain initially, when there is only one point that is refined. However, consistant with previous results, immediately as the point becomes a region of non-zero measure, idealized performance gains drop to $`30\%`$$`10\%`$. Regardless of the refinement, the TR provides a very small performance enhancement compared to that of the spatial refinement.

The reason for the small predicted efficiency gains, consistent with the discussion of the previous section, is that there quickly become more fine blocks than coarse blocks in the simulation. By the last frame shown in Figure 4, there are no blocks being evolved at the the coarsest level of refinement, and indeed 80% of the blocks are at the highest level of refinement. Thus, even if all other blocks required zero work to evolve, we could only achieve a $`20\%`$ speedup.

Next we consider an interface problem – a 2d detonation that will eventually undergo a cellular instability. These simulations are from results published in celldet2d . A mesh is shown in Figure 5. This corresponds almost exactly to the idealized interface problem of the previous section, but here the domain is very long in one direction, increasing the number of low-cost coarsest blocks in the domain. This change in distribution of blocks means that this problem can benefit more from TR. The numerical results are shown in Table 4.

Here we see TR’s efficiency gains actually decrease with increasing resolution, and also see a familiar pattern of TRs efficiency gains going in the opposite direction of spatial AMR efficiency gains. Even at the resolution where TRs efficiency gains are largest, they are much smaller than the improvement from using spatial AMR.

We next consider the development of the Rayleigh Taylor instability. (Figure 6). This is an interface problem, but in this set of simulations, the center region of the box is resolved to ensure resolution of the velocity perturbations in the region near the interface. Because this region is fully refined, many ‘full cost’ finest blocks are added. This decreases the scope of improvement from TR, as seen in Table 5.

## 4 Conclusion

We have considered efficiency gains for time subcycling for explict or local physics. In these cases the work per block is roughly constant. Further, in most cases there are many more fine blocks than coarse blocks — this is due to simple geometry, as a mesh that refines a significant fraction of its domain will be strongly weighted in favour of small blocks, which must be evolved at a small timestep. Thus, Any attempt to improve performance by focusing on the relatively few larger blocks can only reduce a small fraction of the work that needs to be done to evolve the system one timestep. On the other hand, in studies where only a small number of points in a large domain must be fully resolved, there may be significant efficiency gains from TR methods. Some cosmological hydrodynamical simulations normanextreme  are examples of this situation.

We have not considered here accuracy; taking fewer timesteps may increase accuracy with some solvers, although this isn’t clear for moderately time-accurate algorithms having errors of $`O(\mathrm{\Delta }t^p)`$, $`p>1`$; further, the coarsely refined regions which would benefit from the fewer timesteps are presumably coarsely refined because the overall solution quality is less sensitive to the error in those regions than it is to that of the highly refined parts of the domain. We also do not consider global or implicit solves, where the timestepping algorithm in Fig. 1 must be modified. Global or implicit solves will, depending on the methods used, change the amount of work done per block at different levels of refinement, which can change the results given here considerably.

We have modelled only computational cost in this work. Most of the other costs, cf. §2, work to decrease the efficiency gains of TR. One unmodelled effect that could increase the gains is the reduction of guardcell fills on large blocks. For the oct-tree mesh, where the number of zones per block is fixed, the reduction in guardcell filling work is reduced in the same way as the computational work, so that our conclusions are unchanged. For the patch-based mesh, the effect on the guardcell filling will be dependant on the shape of the refined region and the algorithm used for merging patches of the same refinement level, so that it is difficult to say anything in general.

Thus, block-structured TR significantly enhances performance of local or explicit physics solvers only under fairly narrow circumstances. In circumstances where TR is unlikely to produce much performance enhancement, the added code complexity, memory overhead, and parallel load-balancing issues may make the costs of the technique exceed its benefits.

The authors thank B. Fryxell for useful discussions with this paper, and K. Olson with his help with Paramesh over the past years. We thank A. Calder for data from RT simulations, and F. X. Timmes for data from cellular detonation simulations. We thank T. Plewa, G. Weirs, R. Kirby, and R. Loy for suggesting this work. Support for this work was provided by the Scientific Discovery through Advanced Computing (SciDAC) program of the DOE, grant number DE-FC02-01ER41176 to the Supernova Science Center/UCSC. LJD was supported by the Department of Energy Computational Science Graduate Fellowship Program of the Office of Scientific Computing and Office of Defense Programs in the Department of Energy under contract DE-FG02-97ER25308.

The Flash code is freely available at http://flash.uchicago.edu/.