text:
# Statistical Error in Particle Simulations of Hydrodynamic Phenomena

## 1 Introduction

Recently much attention has been focused on the simulation of hydrodynamic problems at small scales using molecular simulation methods such as Molecular Dynamics (MD)  or the direct simulation Monte Carlo (DSMC) . Molecular Dynamics is generally used to simulate liquids while DSMC is a very efficient algorithm for simulating dilute gases. In molecular simulation methods the connection to macroscopic observable fields, such as velocity and temperature, is achieved through averaging appropriate microscopic properties. The simulation results are therefore inherently statistical and statistical errors due to finite sampling need to be fully quantified.

Though confidence intervals may be estimated by measuring the variance of these sampled quantities, this additional computation can be burdensome and thus is often omitted. Furthermore, it would be useful to estimate confidence intervals a priori so that one could predict the computational effort required to achieve a desired level of accuracy. For example, it is well known that obtaining accurate hydrodynamic fields (e.g., velocity profiles) is computationally expensive in low Mach number flows so it is useful to have an estimate of the computational effort required to reach the desired level of accuracy.

In this paper we present expressions for the magnitude of statistical errors due to thermal fluctuations in molecular simulations for the typical observables of interest, such as velocity, density, temperature, and pressure. We also derive expressions for the shear stress and heat flux in the dilute gas limit. Both volume averaging and flux averaging is considered, even though the measurement of the shear stress and heat flux through volume averaging is not exact unless the contribution of the impulsive collisions is accounted for. Although we make use of expressions from equilibrium statistical mechanics, the non-equilibrium modifications to these results are very small, even under extreme conditions . This is verified by the good agreement between our theoretical expressions and the corresponding measurements in our simulations.

In addition to direct measurements of hydrodynamic fields, this analysis will benefit algorithmic applications by providing a framework were statistical fluctuations can be correctly accounted for. One example of such application is the measurement of temperature for temperature-dependent collision rates ; another example is the measurement of velocity and temperature for the purpose of imposing boundary conditions . Additional examples include hybrid methods  where the coupling between continuum and molecular fields requires both averaging of finite numbers of particles in a sequence of molecular realizations as well as the generation of molecular realizations, based on continuum fields, with relatively small numbers of particles (e.g., buffer cells).

In section 2 the theoretical expressions for the statistical error due to thermodynamic fluctuations are derived. These expressions are verified by molecular simulations, as described in section 3. The effect of correlations between samples in dilute gases is briefly discussed in section 4 and concluding remarks appear in section 5.

## 2 Statistical error due to thermal fluctuations

### 2.1 Volume-averaged quantities

We first consider the fluid velocity. In a particle simulation, the flow field is obtained by measuring the instantaneous center of mass velocity, $`𝐮`$, for particles in a statistical cell volume. The statistical mean value of the local fluid velocity, $`𝐮_s`$, is estimated over $`M`$ independent samples. For steady flows, these may be sequential samples taken in time; for transient flows these may be samples from an ensemble of realizations. The average fluid velocity, $`𝐮`$, is defined such that $`𝐮_s𝐮`$ as $`M\mathrm{}`$; for notational convenience we also write $`𝐮=𝐮_0`$. Define $`\delta u_xu_xu_{x0}`$ to be the instantaneous fluctuation in the $`x`$-component of the fluid velocity; note that all three components are equivalent. From equilibrium statistical mechanics ,

$$\delta u_x^2=\frac{kT_0}{mN_0}=\frac{a^2}{\gamma \mathrm{Ac}^2N_0}$$
(1)
where $`N_0`$ is the average number of particles in the statistical cell, $`T_0`$ is the average temperature, $`m`$ is the particle mass, $`k`$ is Boltzmann’s constant, $`a`$ is the sound speed, and $`\gamma =c_P/c_V`$ is the ratio of the specific heats. The acoustic number $`\mathrm{Ac}=a/a^i`$ is the ratio of the fluid’s sound speed to the sound speed of a “reference” ideal gas at the same temperature

$$a^i=\sqrt{\gamma kT/m}$$
(2)
Note that this reference ideal gas has a ratio of specific heats ($`\gamma ^i`$) equal to the original fluid specific heat ratio, that is $`\gamma ^i=\gamma `$ as shown in equation (2).

An alternative construction of (1) is obtained from the equipartition theorem

$$\frac{3}{2}kT_0=\frac{1}{2}m|𝐜𝐜_0|^2=\frac{1}{2}m((c_xc_{x0})^2+(c_yc_{y0})^2+(c_zc_{z0})^2)=\frac{3}{2}m(c_xc_{x0})^2$$
(3)
where $`𝐜`$ is the translational molecular velocity. For a non-equilibrium system, this expression defines $`T_0`$ as the average translational temperature. Note that $`𝐮_0=𝐜_0`$ and

$$|\delta 𝐮|^2=\frac{|𝐜𝐜_0|^2}{N_0}=\frac{3kT_0}{mN_0}.$$
(4)
The above expression also reminds us that the expected error in estimating the magnitude of the fluid velocity is $`\sqrt{3}`$ larger than in estimating a velocity component.

We may define a “signal-to-noise” ratio as the average fluid velocity over its standard deviation; from the above,

$$\frac{|u_{x0}|}{\sqrt{\delta u_x^2}}=\mathrm{Ac}\mathrm{Ma}\sqrt{\gamma N_0}$$
(5)
where $`\mathrm{Ma}=|u_{x0}|/a`$ is the local Mach number based on the velocity component of interest. This result shows that for fixed Mach number, in a dilute gas simulation ($`\mathrm{Ac}=1`$), the statistical error due to thermal fluctuations cannot be ameliorated by reducing the temperature. However, when the Mach number is small enough for compressibility effects to be negligible, favorable relative statistical errors may be obtained by performing simulations at an increased Mach number (to a level where compressibility effects are still negligible).

The one-standard-deviation error bar for the sample estimate $`u_x_s`$ is $`\sigma _u=\sqrt{\delta u_x^2}/\sqrt{M}`$ and the fractional error in the estimate of the fluid velocity is

$$E_u=\frac{\sigma _u}{|u_{x0}|}=\frac{1}{\sqrt{MN_0}}\frac{1}{\mathrm{Ac}\mathrm{Ma}\sqrt{\gamma }},$$
(6)
yielding

$$M=\frac{1}{\gamma \mathrm{Ac}^2N_0\mathrm{Ma}^2E_u^2}.$$
(7)
For example, with $`N_0=100`$ particles in a statistical cell, if a one percent fractional error is desired in a $`\mathrm{Ma}=1`$ flow, about $`M=100`$ independent statistical samples are required (assuming $`\mathrm{Ac}1`$). However, for a $`\mathrm{Ma}=10^2`$ flow, about $`10^6`$ independent samples are needed, which quantifies the empirical observation that the resolution of the flow velocity is computationally expensive for low Mach number flows.

Next we turn our attention to the density. From equilibrium statistical mechanics, the fluctuation in the number of particles in a cell is

$$\delta N^2=N^2\frac{kT_0}{V^2}\left(\frac{V}{P}\right)_T=\kappa _TN_0^2\frac{kT_0}{V}$$
(8)
where $`V`$ is the volume of the statistics cell and $`\kappa _TV^1(V/P)_T`$ is the isothermal compressibility. Note that for a dilute gas $`\kappa _T=1/P`$ so $`\delta N^2=N`$ and, in fact, $`N`$ is Poisson random variable. The fractional error in the estimate of the density is

$$E_\rho =\frac{\sigma _\rho }{\rho _0}=\frac{\sigma _N}{N_0}=\frac{\sqrt{\delta N^2}}{N_0\sqrt{M}}=\frac{\sqrt{\kappa _TkT_0}}{\sqrt{MV}}=\frac{\sqrt{\kappa _T/\kappa _T^i}}{\sqrt{MN_0}}$$
(9)
where $`\kappa _T^i=V/N_0kT_0`$ is the isothermal compressibility of the reference dilute gas ($`\gamma ^i=\gamma `$) at the same density and temperature. Since $`a1/\sqrt{\kappa _T}`$,

$$E_\rho =\frac{1}{\sqrt{MN_0}}\frac{1}{\mathrm{Ac}}$$
(10)
Note that for fixed $`M`$ and $`N_0`$, the error decreases as the compressibility decreases (i.e., as the sound speed increases) since the density fluctuations are smaller.

Let us now consider the measurement of temperature. First we should remark that the measurement of instantaneous temperature is subtle, even in a dilute gas. But given that temperature is measured correctly, equilibrium statistical mechanics gives the variance in the temperature fluctuations to be

$$\delta T^2=\frac{kT_0^2}{c_VN_0}$$
(11)
where $`c_V`$ is the heat capacity per particle at constant volume. The fractional error in the estimate of the temperature is

$$E_T=\frac{\sigma _T}{T_0}=\frac{\sqrt{\delta T^2}}{T_0\sqrt{M}}=\frac{1}{\sqrt{MN_0}}\sqrt{\frac{k}{c_V}}$$
(12)
Because the fluctuations are smaller, the error in the temperature is smaller when the heat capacity is large. Note that the temperature associated with various degrees of freedom (translational, vibrational, rotational) may be separately defined and measured. For example, if we consider only the measurement of the translational temperature, then the appropriate heat capacity is that of an ideal gas with three degrees of freedom, i.e. $`c_V=\frac{3}{2}k`$, corresponding to the three translational components.

Finally, the variance in the pressure fluctuations is

$$\delta P^2=kT_0\left(\frac{P}{V}\right)_S=\frac{\gamma kT_0}{V\kappa _T}$$
(13)
so the fractional error in the estimate of the pressure is

$$E_P=\frac{\sigma _P}{P_0}=\frac{\sqrt{\delta P^2}}{P_0\sqrt{M}}=\frac{P_0^i}{P_0}\frac{\mathrm{Ac}\sqrt{\gamma }}{\sqrt{MN_0}}$$
(14)
where $`P_0^i=N_0kT_0/V`$ is the pressure of an ideal gas under the same conditions. Note that the error in the pressure is proportional to the acoustic number while the error in the density, eqn. (10), goes as $`\mathrm{Ac}^1`$.

### 2.2 Shear stress and heat flux for dilute gases

The thermodynamic results in the previous section are general; in this section we consider transport quantities and restrict our analysis to dilute gases. In a dilute gas, the shear stress and heat flux are defined as

$$\tau _{xy0}=\tau _{xy}=\rho c_xc_y$$
(15)
and

$$q_{x0}=q_x=\frac{1}{2}\rho c_xc^2$$
(16)
respectively. This definition neglects the contribution of impulsive collisions to the transport within a given volume due to the negligible size of particles in the dilute limit. When, however, simulation methods such as DSMC are used to simulate a dilute gas that use a finite molecular size, a small error arises when momentum and energy transport is calculated in a volume averaged manner. In fact, this inconsistency has been addressed at the equation of state level by Alexander et al.  who pointed out that DSMC reproduces an ideal gas equation of state while simulating transport of a hard-sphere gas with molecules of finite size. Thus, under the assumption of a very dilute gas, the fluctuation in the (equilibrium) shear stress and heat flux in a volume $`V`$ containing $`N_0`$ particles can be calculated using the Maxwell-Boltzmann distribution. In equilibrium, the expected values of the above quantities are zero.

Using the definitions of shear stress and heat flux in terms of moments of the velocity distribution, direct calculation of the variance of the $`x`$-$`y`$ component of the stress tensor based on a single-particle distribution function gives

$`\tau _{xy}^2`$ $`=`$ $`(\rho c_xc_y)^2`$ (17)
$`=`$ $`\rho _0^2(c_xc_y)^2`$ (18)
$`=`$ $`{\displaystyle \frac{1}{4}}\rho _0^2\left({\displaystyle \frac{2kT_0}{m}}\right)^2=P_0^2`$ (19)
In obtaining the second equation we assumed $`c_x=c_y=(\rho \rho _0)c_xc_y=0`$. For the $`x`$ component of the heat flux vector, we find

$`q_x^2`$ $`=`$ $`({\displaystyle \frac{1}{2}}\rho c_xc^2)^2`$ (20)
$`=`$ $`\rho _0^2({\displaystyle \frac{1}{2}}c_xc^2)^2`$ (21)
$`=`$ $`{\displaystyle \frac{35}{32}}\rho _0^2\left({\displaystyle \frac{2kT_0}{m}}\right)^3={\displaystyle \frac{35}{8}}c_\mathrm{m}^2P_0^2`$ (22)
where $`c_\mathrm{m}=\sqrt{2kT_0/m}`$ the most probable particle speed and we have assumed $`c_x=c_y=(\rho \rho _0)c^2=0`$. Note that in equilibrium for a cell containing $`N_0`$ particles, the variance of the mean is given by $`\delta \tau _{xy}^2=\tau _{xy}^2/N_0`$ and $`\delta q_x^2=q_x^2/N_0`$ for the shear stress and heat flux respectively.

In order to derive expressions for the relative fluctuations we need expressions for the magnitude of the fluxes. We are only able to provide closed form expressions for the latter in the continuum regime where

$$\tau _{xy}=\mu \left(\frac{u_x}{y}+\frac{u_y}{x}\right)$$
(23)
and

$$q_x=\kappa \frac{T}{x}$$
(24)
where $`\mu `$ is the coefficient of viscosity and $`\kappa `$ is the thermal conductivity. Above Knudsen numbers of $`\mathrm{Kn}0.1`$, it is known that these continuum expressions are only approximate and better results are obtained using more complicated formulations from kinetic theory (e.g., Burnett’s formulation). Here, the Knudsen number is defined as $`\mathrm{Kn}=\lambda /\mathrm{}`$ and the mean free path as

$$\lambda =\frac{8}{5\sqrt{\pi }}\frac{c_\mathrm{m}\mu }{P_0}$$
(25)
Note that this expression for the mean free path simplifies to the hard sphere result when the viscosity is taken to be that of hard spheres.

Using (23) and (24) we find that in *continuum flows*, the relative fluctuations in the shear stress and heat flux are given by

$$E_\tau =\frac{\sqrt{\delta \tau _{xy}^2}}{|\tau _{xy0}|\sqrt{M}}=\frac{16}{5\sqrt{2\pi \gamma }}\frac{1}{\mathrm{Kn}\mathrm{Ma}_{}}\frac{1}{\sqrt{N_0M}}\left|\frac{u_x^{}}{y^{}}+\frac{u_y^{}}{x^{}}\right|^1$$
(26)
and

$$E_q=\frac{\sqrt{\delta q_x^2}}{|q_{x0}|\sqrt{M}}=\frac{8\sqrt{35}}{5\sqrt{2\pi }}\frac{\mathrm{Pr}(\gamma 1)}{\gamma }\frac{T}{\mathrm{\Delta }T}\frac{1}{\mathrm{Kn}}\frac{1}{\sqrt{N_0M}}\left|\frac{T^{}}{x^{}}\right|^1$$
(27)
respectively. Here, stars denote non-dimensional quantities: $`u^{}=u_0/\stackrel{~}{u}`$, $`T^{}=T/\mathrm{\Delta }T`$ and $`x^{}=x/\mathrm{}`$, where $`\stackrel{~}{u}`$, $`\mathrm{\Delta }T`$, and $`\mathrm{}`$ are characteristic velocity, temperature variation and length. The Mach number $`\mathrm{Ma}_{}`$ is defined with respect to the characteristic velocity $`\stackrel{~}{u}`$ rather than the local velocity as in eq. (6).

If viscous heat generation is responsible for the temperature differences characterized by $`\mathrm{\Delta }T`$, then it is possible to express equation (27) in the following form

$$E_q=\frac{\sqrt{\delta q_x^2}}{|q_{x0}|\sqrt{M}}=\frac{8\sqrt{35}}{5\gamma \sqrt{2\pi }}\frac{\mathrm{Br}}{\mathrm{KnMa}_{}^2}\frac{1}{\sqrt{N_0M}}\left|\frac{T^{}}{x^{}}\right|^1$$
(28)
The Brinkman number

$$\mathrm{Br}=\frac{\mu \stackrel{~}{u}^2}{\kappa \mathrm{\Delta }T},$$
(29)
is the relevant non-dimensional group that compares temperature differences due to viscous heat generation to the characteristic temperature differences in the flow. (It follows that if viscous heat generation is responsible for the temperature changes, $`\mathrm{Br}1`$.)

It is very instructive to extend the above analysis to equation (12). If we define the relative error in temperature with respect to the temperature changes rather than the absolute temperature, we obtain

$`E_{\mathrm{\Delta }T}`$ $`=`$ $`{\displaystyle \frac{\sigma _T}{\mathrm{\Delta }T}}={\displaystyle \frac{\sqrt{\delta T^2}}{\mathrm{\Delta }T\sqrt{M}}}={\displaystyle \frac{T_0}{\mathrm{\Delta }T\sqrt{MN_0}}}\sqrt{{\displaystyle \frac{k}{c_V}}}`$ (30)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{Pr}(\gamma 1)}}{\displaystyle \frac{\mathrm{Br}}{\mathrm{Ma}_{}^2}}{\displaystyle \frac{1}{\sqrt{MN_0}}}\sqrt{{\displaystyle \frac{k}{c_V}}}`$ (31)
where again, if viscous heat generation is the only source of heat, $`\mathrm{Br}1`$. The above development shows that resolving the temperature differences or heat flux due to viscous heat generation is very computationally inefficient for low speed flow since for a given expected error $`E_{\mathrm{\Delta }T}`$ we find that the number of samples $`M\mathrm{Ma}_{}^4`$.

Comparison of equations (6) and (26) and equations (12) and (27) reveals that

$$E_\tau \frac{E_u}{\mathrm{Kn}}$$
(32)
and

$$E_q\frac{E_{\mathrm{\Delta }T}}{\mathrm{Kn}}$$
(33)
since the non-dimensional gradients will be of order one. As the above equations were derived for the continuum regime ($`\mathrm{Kn}<0.1`$), it follows that the relative error in these moments is significantly higher. This will also be shown to be the case in the next section when the shear stress and heat flux are evaluated as fluxal (surface) quantities. This has important consequences in hybrid methods , as coupling in terms of state (Dirichlet) conditions is subject to less variation than coupling in terms of flux conditions.

### 2.3 Fluxal quantities

In this section we predict the relative errors in the fluxes of mass and heat as well as the components of the stress tensor, when calculated as fluxes across a reference surface. Our analysis is based on the assumption of an infinite, ideal gas in equilibrium. In this case, it is well known that the number of particles, $`N_i`$, in each infinitessimal cell $`i`$ is an independent Poisson random variable with mean and variance,

$$N_i=\delta N_i^2=n\mathrm{\Delta }x_i\mathrm{\Delta }y_i\mathrm{\Delta }z_i,$$
(34)
where $`n=\rho /m`$ is the particle number density.

In the Appendix, it is shown that, more generally, all number fluctuations are independent and Poisson distributed. In particular, the number of particles, $`N_{ij}^+`$, leaving a cell $`i`$ at position $`(x_i,y_i,z_i)`$, crossing a surface, and arriving in another cell $`j`$ at position $`(x_j,y_j,z_j)`$ after a time $`\mathrm{\Delta }t`$ (see Figure 1) is also an independent Poisson random variable. Its mean and variance are given by

$$N_{ij}^+=(\delta N_{ij}^+)^2=P_{ij}N_i$$
(35)
where

$$P_{ij}=P(x_{ij},y_{ij},z_{ij})\mathrm{\Delta }x_j\mathrm{\Delta }y_j\mathrm{\Delta }z_j,$$
(36)
is the probability that a particle makes the transition, which depends only upon the relative displacements, $`x_{ij}=x_jx_i`$, $`y_{ij}=y_jy_i`$, and $`z_{ij}=z_jz_i`$. These properties hold exactly for an infinite ideal gas in equilibrium, but they are also excellent approximations for most finite, dilute gases, even far from equilibrium, e.g. as shown by the DSMC result in Figure 2.

Without loss of generality we place a large planar surface of area $`A`$ at $`x=0`$. The flux in the positive direction, $`J_h^+`$, of any velocity-dependent quantity, $`h(c_x,c_y,c_z)`$, can be expressed as a sum of independent, random contributions from different cells $`i`$ and $`j`$ on opposite sides of the surface,

$$J_h^+=\frac{1}{A\mathrm{\Delta }t}\underset{x_i<0}{}\underset{x_j>0}{}\underset{y_i}{}\underset{y_j}{}\underset{z_i}{}\underset{z_j}{}\overline{h}(x_{ij},y_{ij},z_{ij})N_{ij}^+$$
(37)
where $`\overline{h}(x_{ij},y_{ij},z_{ij})=h(c_x,c_y,c_z)`$ is the same quantity expressed in terms of relative displacements during the time interval, $`\mathrm{\Delta }t`$. We begin by calculating the mean flux. Taking expectations in Eq. (37), we obtain

$`J_h^+`$ $`=`$ $`{\displaystyle \frac{1}{A\mathrm{\Delta }t}}{\displaystyle \underset{x_i<0}{}}{\displaystyle \underset{x_j>0}{}}{\displaystyle \underset{y_i}{}}{\displaystyle \underset{y_j}{}}{\displaystyle \underset{z_i}{}}{\displaystyle \underset{z_j}{}}\overline{h}(x_{ij},y_{ij},z_{ij})N_iP_{ij}`$ (38)
$``$ $`{\displaystyle \frac{1}{\mathrm{\Delta }t}}{\displaystyle \underset{x_i<0}{}}{\displaystyle \underset{x_j>0}{}}{\displaystyle \underset{y_j}{}}{\displaystyle \underset{z_j}{}}\overline{h}(x_{ij},y_j,z_j)n\mathrm{\Delta }x_iP(x_{ij},y_j,z_j)\mathrm{\Delta }x_j\mathrm{\Delta }y_j\mathrm{\Delta }z_j,`$ (39)
where without loss of generality we have taken $`y_i=z_i=0`$ and also separately performed the sums parallel to the surface,

$$\underset{y_i}{}\underset{z_i}{}\mathrm{\Delta }y_i\mathrm{\Delta }z_i=A,$$

by invoking translational invariance and neglecting any edge effects, since $`c_y\mathrm{\Delta }t,c_z\mathrm{\Delta }t=o(\sqrt{A})`$ in the limit $`\mathrm{\Delta }t0`$ (taken below).

Passing to the continuum limit in Eq. (39) and letting $`A\mathrm{}`$, we arrive at an integral expression for the mean flux through an infinite flat surface,

$`J_h^+`$ $`=`$ $`{\displaystyle \frac{n}{\mathrm{\Delta }t}}{\displaystyle _{\mathrm{}}^0}𝑑x^{}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑y{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z\overline{h}(xx^{},y,z)P(xx^{},y,z)`$ (40)
$`=`$ $`{\displaystyle \frac{n}{\mathrm{\Delta }t}}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle _x^{\mathrm{}}}𝑑s{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑y{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z\overline{h}(s,y,z)P(s,y,z).`$
Since the integrand does not depend on $`x`$, switching the order of integration,

$`J_h^+`$ $`=`$ $`{\displaystyle \frac{n}{\mathrm{\Delta }t}}{\displaystyle _0^{\mathrm{}}}𝑑s{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑y{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z\overline{h}(s,y,z)P(s,y,z){\displaystyle _0^s}𝑑x`$ (41)
$`=`$ $`n{\displaystyle _0^{\mathrm{}}}𝑑s{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑y{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z\overline{h}(s,y,z)P(s,y,z){\displaystyle \frac{s}{\mathrm{\Delta }t}}`$
produces a simple formula for the mean flux as a conditional average over the velocity distribution (with $`c_x>0`$),

$$J_h^+=nc_xh(c_x,c_y,c_z)^+$$
(42)
in the limit $`\mathrm{\Delta }t0`$. For $`h(c_x,c_y,c_z)=m,mc_x,mc_y,mc_z,\frac{1}{2}mc^2`$ in the absence of a mean flow, we obtain,

$`J_m^+`$ $`=`$ $`nc_xm^+={\displaystyle \frac{\rho c_\mathrm{m}}{2\sqrt{\pi }}}`$ (43)
$`J_{xx}^+`$ $`=`$ $`nc_xmc_x^+={\displaystyle \frac{\rho kT_0}{2m}}={\displaystyle \frac{P_0}{2}}`$ (44)
$`J_{xy}^+`$ $`=`$ $`nc_xmc_y^+=0J_{xz}^+=nc_xmc_z^+=0`$ (45)
$`J_e^+`$ $`=`$ $`nc_x\frac{1}{2}mc^2^+={\displaystyle \frac{\rho c_\mathrm{m}^3}{2\sqrt{\pi }}}`$ (46)
which are the well-known results for the expected one-sided fluxes of mass, momentum, and energy, respectively.

Using the same formalism as above, we now derive a general formula for the variance of the flux, $`(\delta J_h^+)^2`$. For simplicity, we assume that the mean flow is zero normal to the reference surface, but our results below are easily extended to the case of a non-zero mean flow through the surface since Eq. (35) still holds in that case (see Figure 3). Taking the variance of Eq. (37) and using the independence of $`\{N_{ij}^+\}`$, we obtain

$`(\delta J_h^+)^2`$ $`=`$ $`{\displaystyle \frac{1}{A^2\mathrm{\Delta }t^2}}{\displaystyle \underset{x_i<0}{}}{\displaystyle \underset{x_j>0}{}}{\displaystyle \underset{y_i}{}}{\displaystyle \underset{y_j}{}}{\displaystyle \underset{z_i}{}}{\displaystyle \underset{z_j}{}}\overline{h}(x_{ij},y_{ij},z_{ij})^2(\delta N_{ij}^+)^2`$ (47)
$`=`$ $`{\displaystyle \frac{1}{A^2\mathrm{\Delta }t^2}}{\displaystyle \underset{x_i<0}{}}{\displaystyle \underset{x_j>0}{}}{\displaystyle \underset{y_i}{}}{\displaystyle \underset{y_j}{}}{\displaystyle \underset{z_i}{}}{\displaystyle \underset{z_j}{}}\overline{h}(x_{ij},y_{ij},z_{ij})^2N_iP_{ij}`$
Following the same steps above leading from Eq. (38) to Eq. (42) we arrive at the simple formula,

$$(\delta J_h^+)^2=\frac{n}{A\mathrm{\Delta }t}c_xh(c_x,c_y,c_z)^2^+$$
(48)
in the continuum limit as $`\mathrm{\Delta }t0`$. The variance of the total flux is

$$\delta J_h^2=(\delta J_h^++\delta J_h^{})^2=(\delta J_h^+)^2+(\delta J_h^{})^2=2(\delta J_h^+)^2$$
(49)
since the one-sided fluxes through the surface, $`J_h^+`$ and $`J_h^{}`$, are independent and identically distributed (with opposite sign).

Using this general result, we can evaluate the standard deviations of the fluxes above,

$`\sqrt{\delta J_m^2}`$ $`=`$ $`\sqrt{{\displaystyle \frac{\rho c_\mathrm{m}}{\sqrt{\pi }}}}\sqrt{{\displaystyle \frac{m}{A\mathrm{\Delta }t}}}`$ (50)
$`\sqrt{\delta J_{xy}^2}`$ $`=`$ $`\sqrt{{\displaystyle \frac{\rho c_\mathrm{m}^3}{2\sqrt{\pi }}}}\sqrt{{\displaystyle \frac{m}{A\mathrm{\Delta }t}}}`$ (51)
$`\sqrt{\delta J_e^2}`$ $`=`$ $`\sqrt{{\displaystyle \frac{3\rho c_\mathrm{m}^5}{2\sqrt{\pi }}}}\sqrt{{\displaystyle \frac{m}{A\mathrm{\Delta }t}}}`$ (52)
These formulae may be simplified by noting that

$$J_m^+=\frac{mN^+}{A\mathrm{\Delta }t}=\frac{\rho c_\mathrm{m}}{2\sqrt{\pi }}$$
(53)
where $`N^+`$ is the mean total number of particles crossing the reference surface in one direction in time $`\mathrm{\Delta }t,`$, which yields

$`\sqrt{\delta J_m^2}`$ $`=`$ $`{\displaystyle \frac{\rho c_\mathrm{m}}{\sqrt{2\pi }}}{\displaystyle \frac{1}{\sqrt{N^+}}}`$ (54)
$`\sqrt{\delta J_{xy}^2}`$ $`=`$ $`{\displaystyle \frac{\rho c_\mathrm{m}^2}{2\sqrt{\pi }}}{\displaystyle \frac{1}{\sqrt{N^+}}}`$ (55)
$`\sqrt{\delta J_e^2}`$ $`=`$ $`\sqrt{{\displaystyle \frac{3}{\pi }}}{\displaystyle \frac{\rho c_\mathrm{m}^3}{2}}{\displaystyle \frac{1}{\sqrt{N^+}}}`$ (56)
We may relate these to the results from the previous section by identifying $`\delta (\tau _{xy}^f)^2=\delta J_{xy}^2`$ and $`\delta (q_x^f)^2=\delta J_e^2`$, and additionally include the effect of $`M`$ (independent) samples in time. The superscript $`f`$ denotes fluxal measurement. By noting that transport fluxes are defined with respect to the rest frame of the fluid, it can be easily verified that the above relations hold in the case where a mean flow in directions parallel to the measuring surface exists, under the assumption of a local equilibrium distribution.

We can derive expressions for the relative expected error in the continuum regime in which models exist for the shear stress and heat flux. In this regime we find

$`E_\tau ^f={\displaystyle \frac{\sqrt{\delta (\tau _{xy}^f)^2}}{|\tau _{xy0}^f|\sqrt{M}}}`$ $`=`$ $`{\displaystyle \frac{16}{5\pi \sqrt{2\gamma }}}{\displaystyle \frac{1}{\mathrm{Kn}\mathrm{Ma}_{}}}{\displaystyle \frac{1}{\sqrt{MN^+}}}\left|{\displaystyle \frac{u^{}}{y^{}}}+{\displaystyle \frac{v^{}}{x^{}}}\right|^1`$ (57)
and

$`E_q^f={\displaystyle \frac{\sqrt{\delta (q_x^f)^2}}{|q_{x0}^f|\sqrt{M}}}`$ $`=`$ $`{\displaystyle \frac{16\sqrt{3}}{5\pi }}{\displaystyle \frac{\mathrm{Pr}(\gamma 1)}{\gamma }}{\displaystyle \frac{T}{\mathrm{\Delta }T}}{\displaystyle \frac{1}{\mathrm{Kn}}}{\displaystyle \frac{1}{\sqrt{MN^+}}}\left|{\displaystyle \frac{T}{x^{}}}\right|^1`$ (58)
$`=`$ $`{\displaystyle \frac{16\sqrt{3}}{5\gamma \pi }}{\displaystyle \frac{\mathrm{Br}}{\mathrm{Kn}\mathrm{Ma}_{}^2}}{\displaystyle \frac{1}{\sqrt{MN^+}}}\left|{\displaystyle \frac{T}{x^{}}}\right|^1`$ (59)
Comparing (57) with the corresponding expressions for volume-averaged stress tensor, (28), one finds that, aside from the numerical coefficients, the expressions differ only in the number of particles used, either $`N^+`$ or $`N_0`$; one finds a similar result for the heat flux.

### 2.4 Connection to Fluctuating Hydrodynamics

Fluctuating hydrodynamics, as developed by Landau, approximates the stress tensor and heat flux as white noises, with variances fixed by matching equilibrium fluctuations. In this section we identify the connection between Landau’s theory and the variances of fluxes obtained in section 2.2.

Landau introduced fluctuations into the hydrodynamic equations by adding white noise terms to the stress tensor and heat flux  (in the spirit of Langevin’s theory of Brownian motion ). The amplitudes of these noises are fixed by evaluating the resulting variances of velocity and temperature and matching with the results from equilibrium statistical mechanics . For example, in Landau’s formulation the total heat flux in the $`x`$-direction is $`q_x^\mathrm{L}=\kappa T/x+g_x`$ where the first term on the r.h.s. is the deterministic part of the flux and the second is the white noise term. The latter has mean zero and time correlation given by the following expression

$$g_x(t)g_x(t^{})=\frac{2k\kappa T^2}{V}\delta (tt^{})$$
(60)
Note that at a steady state the deterministic part is constant so $`\delta q_x^\mathrm{L}(t)\delta q_x^\mathrm{L}(t^{})=g_x(t)g_x(t^{})`$.

On the other hand, recall from eqn. (22),

$$\delta q_x^2=\frac{35}{32}\frac{\rho ^2c_\mathrm{m}^6}{N_0}$$
(61)
The question naturally arises: how does one reconcile (60) and (61)? Note that the fluctuating hydrodynamics expression contains the thermal conductivity, which depends on the particle interaction (e.g., for hard-spheres $`\kappa `$ depends on the particle diameter) while the kinetic theory expression is independent of this interaction.

The key lies in identifying the $`\delta `$-function with a decay time, $`t_d`$, that is vanishingly small at hydrodynamic scales. Specifically, we may write

$$\delta (tt^{})\{\begin{array}{cc}t_d^1& |tt^{}|<t_d\\ 0& \mathrm{otherwise}\end{array}$$
(62)
so

$$g_x^2(t)=\frac{2k\kappa T^2}{t_dV}=\frac{\kappa V}{2kN_0c_\mathrm{m}^2t_d}\frac{\rho ^2c_\mathrm{m}^6}{N_0}$$
(63)
Comparing with the above gives

$$t_d=\frac{16}{35}\frac{\kappa V}{kN_0c_\mathrm{m}^2}$$
(64)
For a hard sphere gas, $`\kappa =\frac{5}{2}c_V\mu `$ so using (25) we may write this as,

$$t_d=\frac{15\sqrt{\pi }}{28}\frac{\lambda }{c_\mathrm{m}}$$
(65)
For other particle interactions the coefficients will be slightly different but in general $`t_d\lambda /c_\mathrm{m}`$, thus it is approximately equal to the molecular collision time. In conclusion, the two formulations are compatible once the white noise approximation in fluctuating hydrodynamics is justified by the fact that the hydrodynamic time scale is much longer than the kinetic (i.e., collisional) time scale. Landau’s construction provides a useful hydrodynamic approximation for $`g_x`$ but (61) is the actual variance of the heat flux.

## 3 Simulations

### 3.1 Dilute Gases

We performed DSMC simulations to verify the validity of the expressions derived above. Standard DSMC techniques  were used to simulate flow of gaseous argon (molecular mass $`m=6.63\times 10^{26}`$ kg, hard sphere diameter $`\sigma =3.66\times 10^{10}`$ m) in a two-dimensional channel (length $`L`$ and height $`H`$). The simulation was periodic in the $`x`$ direction (along the channel axis). The two walls at $`y=H/2`$ and $`y=H/2`$ were fully accommodating and flat. The simulation was also periodic in the third (homogeneous) direction.

The average gas density was $`\rho _0=1.78\mathrm{kg}/\mathrm{m}^3`$ and in all calculations over 40 particles per cell were used. The cell size was $`\mathrm{\Delta }x=\lambda _0/3`$ where $`\lambda _0`$ is the reference mean free path. The time step was $`\mathrm{\Delta }t=\lambda _0/(7.5c_\mathrm{m})`$. For a discussion of the errors resulting from finite cell sizes and time steps see . The fractional error in the simulations is obtained from the standard deviation of cell values in the $`x`$ and $`z`$ directions. To ensure that the samples were independent, samples were taken only once every 250 time steps. To ensure that the system was in its steady state the simulation was run for $`10^6`$ time steps before sampling was started.

A constant acceleration was applied to the particles to produce Poiseuille flow in the $`x`$ direction with maximum velocity at the centerline $`u_0^{\mathrm{max}}2`$ m/s. Figures 4, 5, and 6 show good agreement between the theoretical expressions from section 2 and simulation measurements for the fractional error in velocity, density and temperature, respectively. The fractional error in the velocity measurement is minimum at the centerline since the Poiseuille velocity profile is parabolic and maximum at the centerline, (see Fig. 4). The density and temperature were nearly constant across the system so the fractional errors in these quantities are also nearly constant.

The expressions for shear stress and heat flux were verified using Couette (walls at equal temperature with different velocities) and “temperature” Couette (walls at zero velocity with different temperatures) calculations respectively. In these calculations, very small cell sizes ($`\mathrm{\Delta }x=\lambda /6`$) and time steps $`(\mathrm{\Delta }t=\lambda /(30c_\mathrm{m}))`$ were used in order to minimize the discrepancy between the volume-averaged and surface-averaged shear stress and heat flux . The system was equilibrated for $`10^6`$ time steps and samples were taken every 50 time steps. The momentum and energy fluxes de-correlate faster than the conserved hydrodynamic variables, such as density, so independent samples are obtained after fewer time steps (see section 4). Good agreement is found between the theoretical results and simulation measurements for volume averaged and fluxal quantities, as shown in figures 7, 8, 9, 10.

A final note: In DSMC simulations one considers each particle as “representing” a large number of molecules in the physical system. In all the expressions given above, $`N_0`$ and $`N^+`$ relates to the number of particles used by the simulation so the fluctuations can be reduced by using larger numbers of particles (i.e., using a lower molecule-to-particle ratio).

### 3.2 Dense fluids

We performed molecular dynamics simulations to test the validity of equations (6), (10), (12) for dense fluids. A similar geometry to the dilute gas simulations described above was used but at a significantly higher density. In particular, we simulated liquid argon ($`\sigma _{LJ}=3.4\times 10^{10}`$m, $`\epsilon _{LJ}=119.8k_b`$) at $`T=240K`$ and $`\rho =860Kg/\mathrm{m}^3`$ in a two-dimensional channel with the $`x`$ and $`z`$ directions periodic. The channel height was $`H=69.7\sigma _{LJ}`$. The wall molecules were connected to fcc lattice sites through springs and interacted with the fluid through a Lennard-Jones potential with the same parameters. The spring constant $`k_s=460\epsilon \sigma ^2`$ was chosen in the way that root mean square displacement of wall atoms around their equilibrium position at the simulated temperature was well below the Lindermann criterion for the melting point of a solid. The length and depth of the system was $`28\sigma _{LJ}`$ and $`29.1\sigma _{LJ}`$ in the $`x`$ and $`z`$ directions respectively. A constant force $`f=8\times 10^5\epsilon /\sigma _{LJ}`$ per particle was used to generate a velocity field with a maximum velocity of approximately $`13`$ m/s.

In order to calculate the fluctuation of density, temperature and velocity, we divided the simulation cell into 13 layers in the $`y`$direction with a height $`\mathrm{\Delta }y=4.6396\sigma `$. We further divided each layer into $`49`$ cells, 7 in each of the $`x`$ and $`z`$ directions. The density, temperature and velocity in each cell were calculated in every $`2000(0.005t_{LJ})`$, where $`t_{LJ}=\sqrt{m_{LJ}\sigma _{LJ}^2/\epsilon _{LJ}}`$. We have checked that this time interval is longer than the system’s correlation time such that samples taken between such intervals are independent. For each cell, $`200`$ samples are used to calculate the average density, temperature and velocity. The fluctuation was calculated for each layer using the 49 equivalent cells in the $`xz`$ plane.

Due to the sensitivity of the compressibility $`\kappa _T`$ on the interaction cutoff $`r_c`$, a rather conservative value of $`r_c=4.0\sigma _{LJ}`$ was used. We also introduced a correction for the still-finite cutoff which used the compressibility predictions of the Modified Benedict-Webb-Rubin equation of state . The agreement between the theoretical predictions and the simulations is good (see Figures 11, 12 and 13).

## 4 Independent Samples and Correlations

The results in figures 7, 8, 9, 10 suggest that volume-averaged measurements provide a superior performance due to a smaller relative error. This conclusion, however, is not necessarily correct because our results are based on arbitrary choices of “measurement spacing”, in the sense that the only consideration was to eliminate correlations in the data since the theoretical formulation in section 2 is based on the assumption of uncorrelated samples. The two methods of sampling are in fact linked by a very interesting interplay between the roles of time and space: fluxal sampling is a measurement at a fixed position in space for a period of time, whereas volume sampling is performed over some region of space at a fixed time. The theoretical performance of each method can be increased by extending the respective window of observation. However, by increasing the time of observation, a fluxal measurement becomes correlated with neighboring measurements if the same particle crosses more than one measuring station within the period of observation. Similarly, by increasing the region of measurement, subsequent volume measurements will suffer from time correlations if previously interrogated particles do not have sufficient time to leave the measurement volume. This relation between spatial and temporal sampling and the role of the particle characteristic velocity is also manifested in the theoretical predictions. Enforcing equality of the variances of the respective volume and fluxal measurements (eqs (55) and (26), and eqs (56) and (27)) yields $`\mathrm{\Delta }x=\beta c_\mathrm{m}\mathrm{\Delta }t`$ where $`\beta 1`$. The generalization of this work to include time and spatial correlations is the subject of future work.

The effect of time correlation in volume measurements can be approximated using the theory of the “persistent random walk” , first introduced by Fürth  and Taylor . A persistent random walk is one in which each step displacement, $`u(t_\eta )=u(\eta \mathrm{\Delta }t)`$, is identically distributed and has a positive correlation coefficient with the previous step,

$$\alpha =\frac{\delta u(t_\eta )\delta u(t_{\eta 1})}{\delta u^2},$$
(66)
($`0<\alpha <1`$), e.g. to model diffusion in a turbulent fluid . This assumption implies that step correlations decay exponentially in time,

$$\frac{\delta u(t_\eta )\delta u(0)}{\delta u^2}=\alpha ^\eta =e^{t_\eta /t_c},$$
(67)
where $`t_c=\mathrm{\Delta }t/\mathrm{log}\alpha `$ is the correlation time, beyond which the steps are essentially indepedent. The position, $`U(t_\eta )`$, of the random walker after $`\eta `$ steps is the sum of these correlated random displacements, $`U(t_\eta )=_{j=1}^\eta u(t_j)`$. Following Taylor, it is straightforward to show that for times long compared to the correlation time ($`t_\eta t_c`$), the usual diffusive scaling holds

$$\delta U(t_\eta )^2\delta u^2\eta \left(\frac{1+\alpha }{1\alpha }\right),$$
(68)
where the bare diffusion coefficient is modified by the term in parentheses due to correlations.

There is a natural connection with the measurement of statistical averages. If we view each step in the persistent random walk as a correlated sample of some quantity in the gas, then the position of the walker (divided by the number of samples) corresponds to the sample average. Thus the variance of a set of sequentially correlated random variables $`\{u(t_\eta )\}`$, $`Var^c(u)`$, may be written as

$$Var^c(u)=Var(u)\frac{1+\alpha }{1\alpha }$$
(69)
where $`Var(u)`$ is the variance of the uncorrelated samples ($`\alpha =0`$). The theory above implies that the sample variance is amplified by the presence of correlations, because effectively fewer independent samples have been taken compared to the uncorrelated case, $`\alpha =0`$.

Note that a sequence of correlated random variables, $`\{u(t_\eta )\}`$, satisfying Eqs. (66) and (67) can be explicitly constructed from a sequence of independent, identically distributed variables, $`\{\stackrel{~}{u}(t_\eta )\}`$, by letting $`u(t_\eta )`$ equal the previous value $`u(t_{\eta 1})`$ with probability $`\alpha `$ or a new value $`\stackrel{~}{u}(t_\eta )`$ with probability $`1\alpha `$. This allows us to interpret $`\alpha `$ as the probability that a sample is the same as the previous one, which is precisely the source of correlations when volume-sampling dilute gases.

There are two distinct ways that a new sample of a dilute gas can actually provide new information: (i) Either some new particles have entered the sampling cell, or (ii) some particles previously inside the cell have changed their properties due to collisions. Regarding (i), the probability that a particle will remain in a cell of size $`\mathrm{\Delta }x`$ after time step $`\mathrm{\Delta }t`$ is given by

$$\alpha _l=\left(\frac{1}{\mathrm{\Delta }x}_{\mathrm{\Delta }x/2}^{\mathrm{\Delta }x/2}_{\mathrm{\Delta }x/2}^{\mathrm{\Delta }x/2}f\left(\frac{x\stackrel{~}{x}}{\mathrm{\Delta }t}\right)𝑑x𝑑\stackrel{~}{x}\right)^d$$
(70)
where $`d`$ is the dimensionality of the cell and $`f(c_i)`$ is the probability distribution function of the particle velocity, and thus $`f(\frac{x\stackrel{~}{x}}{\mathrm{\Delta }t})`$ represents the probability that a particle originating from location $`\stackrel{~}{x}`$ is found at $`x`$ after a time interval, $`\mathrm{\Delta }t`$. Note that the above expression holds when the components of the particle velocity in different directions are uncorrelated and that the particle spatial distribution inside the cell is initially uniform. Regarding (ii), if the single-particle auto-correlation function decays exponentially, the “renewal” probability $`1\alpha _c`$ for collisional de-correlation can, at least in principle, be inferred from simulations using Eq. (67). The net correlation coefficient, giving the probability of an “identical” sample, is $`\alpha =\alpha _c\alpha _l`$.

In our DSMC simulations, the single-particle correlation coefficient for the velocity was estimated to be $`\alpha _c0.94`$ for $`\mathrm{\Delta }t=\lambda /(15c_\mathrm{m})`$ by fitting Eq. (67) to data; the same value for $`\alpha _c`$ was found in both equilibrium and non-equilibrium simulations. In the comparison of figures 14, 15 and 16, however, we use $`\alpha _c=1`$ since mass, momentum and energy are conserved during collisions. We find that this analysis produces very good results in the case of density and temperature (see figures 15 and 16) and acceptable results for the case of mean velocity (see figure 14). Use of $`\alpha _c=0.94`$ would tend to make the agreement better but no theoretical justification exists for it. The value of $`\alpha _l`$ was directly calculated by assuming an equilibrium distribution. For our two-dimensional calculations with $`\mathrm{\Delta }x/(c_\mathrm{m}\mathrm{\Delta }t)=2.5`$, we find $`\alpha _l0.6`$. The effect of a mean flow is very small if the Mach number is small. This was verified through both direct evaluation of Eq. (70) using a local equilibrium distribution function and DSMC simulations.

## 5 Conclusions

We have presented expressions for the statistical error in estimating the velocity, density, temperature and pressure in molecular simulations. These expressions were validated for flow of a dilute gas and dense liquid in a two-dimensional channel using the direct simulation Monte Carlo and Molecular Dynamics respectively. Despite the non-equilibrium nature of the validation experiments, good agreement is found between theory and simulation, verifying that modifications to non-equilibrium results are very small. In particular, in the dense fluid case, despite the significant non-equilibrium due to a shear of the order of $`5.5\times 10^8\mathrm{s}^1`$, the agreement with equilibrium theory is remarkable. We thus expect these results to hold for general non-equilibrium applications of interest.

Predictions were also presented for the statistical error in estimating the shear stress and heat flux in dilute gases through cell averaging and surface averaging. Comparison with direct Monte Carlo simulations shows that the equilibrium assumption is justified. Within the same sets of assumptions, we were able to show that the distribution of particles leaving a cell is Poisson-distributed.

One consideration that significantly limits the applicability of cell (volume) averaging for transport quantities is their neglect of transport due to collisions in the cell volume. In DSMC in particular, as the time step of the simulation is increased, the fraction of particles in a cell undergoing collision increases, and as a result, the error from using the above method increases.

It was found that the fluctuation in state variables is significantly smaller compared to flux variables in the continuum regime $`Kn0`$. This is important for the development of hybrid methods. Although a direct comparison was only presented between volume averaged quantities, we find that the fluxal measurements for the shear stress and heat flux perform similarly to the volume averaged counterparts (regarding scaling with the Knudsen number).

## Appendix: Number Fluctuations in Dilute Gases

Consider an infinite ideal gas in equilibrium with mean number density, $`n`$. By definition, each infinitessimal volume element, $`dV`$, contains a particle with probability, $`ndV`$, and each such event in independent. From these assumptions, it is straightforward to show that the number of particles $`N`$ in an arbitrary volume $`V`$ is a Poisson random variable ,

$$\text{Prob}(N=m)=\frac{e^NN^m}{m!}$$

with mean and variance given by

$$N=\delta N^2=nV.$$

Note that this result does not strictly apply to a finite ideal gas because the probabilities of finding particles in different infinitessimal volumes are no longer independent, due to the global constraint of a fixed total number of particles. Nevertheless, it is an excellent approximation for most dilute gases, even finite non-ideal gases far from equilibrium (e.g. as demonstrated by simulations in the main text).

Now consider a property $`A`$ which each particle in a volume $`V`$ may possess independently with probability, $`P_A`$. In this Appendix, we show that the distribution of the number of such particles, $`N_A`$, is also a Poisson random variable with mean and variance given by

$$N_A=\delta N_A^2=P_AN=P_AnV.$$

For example, in the main text we require the number of particles, $`N_{ij}`$, which travel from one region, $`i`$, to another region, $`j`$, in a time interval, $`dt`$. The proof given here, however, is much more general and applies to arbitrary number fluctuations of an infinite ideal gas in equilibrium, such as the number of particles in a certain region, moving in a certain direction, of a certain “color”, with speeds above a certain theshold, etc.

We begin by expressing $`N_A`$ as a random sum of random variables,

$$N_A=\underset{i=1}{\overset{N}{}}\chi _i$$
(71)
where

$$\chi _i=\{\begin{array}{cc}\hfill 1& \mathrm{if}A\mathrm{occurs}\hfill \\ \hfill 0& \mathrm{otherwise}.\hfill \end{array}$$

is an indicator function for particle $`i`$ to possess property $`A`$, which is a Bernoulli random variable with mean, $`P_A`$. It is convenient to introduce probability generating functions,

$$f_\chi (z)=\underset{m=0}{\overset{\mathrm{}}{}}\text{Prob}(\chi =m)z^m=(1P_A)+P_Az$$

and

$$f_N(z)=\underset{m=0}{\overset{\mathrm{}}{}}\text{Prob}(N=m)z^m=\underset{m=0}{\overset{\mathrm{}}{}}\frac{e^NN^mz^m}{m!}=e^{N(z1)}$$
(72)
because the generating function for a random sum of random variables, as in Eq. (71), is simply given by a composition of the generating functions for the summand and the number of terms ,

$$f_{N_A}(z)=\underset{m=0}{\overset{\mathrm{}}{}}\text{Prob}(N_A=m)z^m=f_N\left(f_\chi (z)\right).$$

Combining these expressions we have

$$f_{N_A}(z)=e^{N(1P_A+P_Az1)}=e^{P_AN(z1)}.$$

Comparing with Eq. (72) completes the proof that $`N_A`$ is a Poisson random variable with mean, $`P_AN`$.

## 6 Acknowledgements

The authors wish to thank M. Malek-Mansour and B. Alder for helpful discussions. The authors would also like to thank X. Garaizar for making this work possible through the computer resources made available to them. This work was also supported, in part, by a grant from the University of Singapore, through the Singapore-MIT alliance.


format: text

text:
# Reflection at large distance from the central engine in Seyferts

## 1 Introduction

In radio-quiet active galactic nuclei, large amounts of cold material are thought to reside at large distance from the central engine. In particular, the unified scheme (Antonucci 1993) postulates a large scale torus (with radius larger than $`10^{16}`$ cm). Such distant cold material (hereafter DCM) may imprint strong reflection features in the hard X-ray spectrum of Seyfert galaxies. These spectral characteristics take the form of a Compton hump in the X-ray spectrum (e.g. Lightman & White 1988; George & Fabian 1991) and a narrow iron line component at 6.4 keV. They add to similar features likely to form much closer to the central source presumably by reflection of the primary X-ray radiation on the inner parts of the accretion disc.

In Compton thick Seyfert 2s, the observed reflection dominated spectra are generally interpreted by DCM reflection with the primary emission obscured presumably by the same outer material (e.g. Reynolds et al. 1994; Matt et al. 2000). In Compton thin Seyfert 2s, Risaliti (2002) recently present, from the analysis of a large number of BeppoSAX observations, strong evidences supporting the presence of a distant reflector in these objects. The unifications models postulate that such material is also present in Seyfert 1s but out of the line of sight. In this context, it has been argued that the reflection produced by an irradiated torus would be sufficient to account for the typical observed reflection spectra in Seyfert 1s (Ghisellini, Haardt & Matt 1994; Krolik, Madau & Życki 1994), without need for reflecting material in the inner parts of the accretion flow.

The first evidence for remote reflection in a type 1 Seyfert came from the BeppoSAX observation of NGC 4051 in a ultra faint state where the nuclear emission had switched-off leaving only a spectrum of pure reflection from material at distance larger than 1017 cm (Guainazzi et al. 1998; Uttley et al. 1999). However, the presence of a reflecting torus in all Seyfert 1s is still a matter of debate.

This important issue motivated the search for a narrow iron emission line that would be the signature for reflection on distant cold and Compton-thick material. Recent Chandra (e.g. Kaspi et al. 2001; Yaqoob et al. 2001) and XMM-Newton (e.g. Reeves et al. 2001; Pounds et al. 2001; Gondoin et al. 2001a, 2001b; Matt et al. 2001, O’Brien et al. 2001) observations indicate that a narrow line is often, if not always, present in Seyfert 1 galaxies, alone or together with a relativistic line component (Tanaka et al. 1995; Fabian et al 2000 and reference therein). An analysis of the composite ASCA spectrum of Seyfert 1s (Lubiñski & Zdziarski 2001, hereafter LZ01) also suggests the presence of a narrow component confirming previous findings on individual sources.

An important complementary approach in order to locate the place of emission of the different components is the study of the variability of the reflected features. Indeed, the temporal response of the line to fluctuations of the continuum depends on the size and distance of the reflector. Reflection features formed in the disc, in the broad line region, or in a torus are expected to show little variability on time scales shorter than a few hours, a few weeks or a few months respectively. In this framework the monitoring campaigns with RXTE of several objects gave indications for a constant line and reflection component flux while the continuum was strongly variable (Papadakis et al. 2002).

There is thus growing evidence in favor of a significant fraction of the reflection features being formed at large distance from the central engine.

On the other hand there is an important observational and theoretical issue that is the existence of a correlation between the amplitude of the reflection component $`R`$ and the spectral index $`\mathrm{\Gamma }`$ found by Zdziarski, Lubiñski & Smith 1999 (hereafter ZLS99) in a large sample of Ginga spectra. The measured $`R`$ tends to be larger in softer sources. This correlation is observed in sample of sources as well as in the time evolution of individual sources. The $`R`$-$`\mathrm{\Gamma }`$ correlation is also found in the BeppoSAX data (Matt 2001), the BeppoSAX sample showing however a significantly different shape, with on average higher $`R`$ values at low $`\mathrm{\Gamma }`$ and a flatter correlation. The usual interpretation of the correlation invokes the feedback from reprocessed radiation emitted by the reflector itself. It thus requires the main source of reflection to be present in the direct environement of the $`X`$-ray source. This thus appears inconsistent with the DCM picture.

One of the main goal of this paper is to attempt solving this paradox. We will show that under assumptions based on the present knowledge of the Seyfert 1s phenomenology, a $`R`$-$`\mathrm{\Gamma }`$ correlation from a distant reflector is actually expected. Another possibility is that the correlation could be spurious. The $`R`$-$`\mathrm{\Gamma }`$ correlation has been indeed widely controversed, in particular regarding the uncertainties on the measurement of $`R`$, but no-one could prove unambiguously that it is an artifact. In this paper we demonstrate that the determination of $`R`$ is sensitive to the assumed high energy cut-off and this may strongly alter the shape of the correlation.

In section 2, we recall the generic effects of the DCM reflection on the measured amplitude of the reflection features. In section 3, we investigate quantitatively the dependence of the reflection amplitude and the iron line equivalent width (EW) on the spectral index $`\mathrm{\Gamma }`$ in the case where the variability of the primary component can be described as a pivoting power-law. Finally, in section 4 we attempt to test the presence of DCM reflection in NGC 5548 using RXTE data.

## 2 Qualitative effects of a distant reflector

As the DCM line is spectrally distinct from the broader line emitted close to the black hole, the existence and the distance/scale of the DCM line region can, in principle, be probed using the new high resolution X-ray telescopes, as discussed above. On the other hand, the DCM Compton hump is more difficult to measure accurately and disentangle from its similar disc counterpart.

Actually, when analyzing the X-ray data, the possible contribution from DCM reflection is usually not considered a priori. The data are generally interpreted in the framework of reflection in the vicinity of the hard X-ray source. Usually, in spectral fits, the shape of the reflection component is computed assuming that the *observed* primary spectrum illuminates an infinite slab (e.g. pexrav model in xspec, Magdziarz & Zdziarski 1995). The normalization of the reflected spectrum is then tuned in order to fit the observed spectrum. The result of this fitting procedure provides the reflection amplitude $`R`$. $`R`$ is normalized so that $`R=1`$ in the case of an isotropic source above an infinite reflecting plane. $`R`$ is often considered as an estimate of $`\mathrm{\Omega }/2\pi `$ where $`\mathrm{\Omega }`$ is the solid angle subtended by the reflector as seen from the isotropic X-ray source. Obviously, any contribution from a remote structure leads to an increase of $`R`$ and may lead to an overestimate of the disc reflection. This may explain the very large $`R`$ coefficients $`2`$, measured in some Seyfert 1s, which are difficult to reconcile with disc reflection (see however Beloborodov 1999; Malzac, Beloborodov & Poutanen 2001).

The nuclei of Seyfert galaxies are known to harbor a significant flux and spectral variability on very short time-scales (see e.g. the recent review by Nandra 2001). Due to its extended structure, the remote reflector cannot respond to the rapid fluctuations of the primary X-ray flux. *The reflected component from the DCM is thus likely to correspond to the time-averaged incident flux seen by the cold material rather than to the instantaneous (i.e. observed) one.* Thus flux changes may induce a significant variation in the $`R`$ value derived from the spectral fits. A flux lower than the average enhances the apparent reflection, on the other hand, a larger flux may reduce $`R`$ down to zero. Then values of $`R`$ as low or large as required by the data can be easily produced. Similarly, the DCM iron line flux is expected to be constant and, consequently, its equivalent width would vary like the inverse of the primary continuum flux around the line energy.

As quoted in the introduction, a significant contribution from a remote reflector seems, at first sight, in contradiction with the reported correlation between $`R`$ and the spectral slope $`\mathrm{\Gamma }`$. Indeed, in the context of DCM reflection, it is difficult to understand why the reflection contribution from DCM should be more important in objects with softer spectra. However, beside the simple effects of flux variability at constant spectrum discussed above, the presence of *spectral* variability also alters the appearance of the DCM reflection. Spectral variations with $`\mathrm{\Delta }\mathrm{\Gamma }`$$`>`$0.1 are indeed commonly reported in Seyfert 1s. In the case of reflection on DCM, such spectral changes have a significant impact on the measured amplitude of the reflection features and the measurement of $`R`$. The resulting $`R`$ vs $`\mathrm{\Gamma }`$ relation then depends on the specific spectral variability mode of the sources that may, under certain circumstances, produce a $`R`$-$`\mathrm{\Gamma }`$ correlation (Nandra et al. 2000).

It appears indeed that the observed variability in Seyfert galaxies is generally consistent with fluctuations of the X-ray spectral slope with a nearly constant comptonised luminosity, i.e. the spectrum is mainly pivoting. This can happen for example if the UV luminosity entering the hot Comptonising plasma changes with a constant heating rate in the hot plasma (e.g. Malzac & Jourdain 2000). This kind of behaviour is widely observed in numerous sources such as NGC 5548 (Nicastro et al. 2000; Petrucci et al. 2000), 3C 120 (Zdziarski & Grandi 2001), NGC 7469 (Nandra et al. 2000; Nandra & Papadakis 2001), 1H 0419-577 (Page et al. 2002), NGC 3783 (De Rosa et al. 2002). In general, the pivot point appears to be located in the 2-10 keV range (see fig. 9 of Petrucci et al. 2000, fig. 3 of Page et al. 2002, fig. 3 of Zdziarski & Grandi 2001).

In the case of a pivoting spectrum at energies lower than 10 keV – i.e. below the energy range were most of reflection is produced – when the spectrum is hard the observed primary flux in the 10–30 keV band is larger, and the relative amplitude of the DCM reflection, $`R`$, is reduced. On the other hand when the spectrum is soft the primary flux in this energy range is enhanced and the measured $`R`$ is larger. This thus produces *a positive correlation between $`R`$ and $`\mathrm{\Gamma }`$*. The lower is the energy of the pivoting point, the steeper the correlation will be.

## 3 Quantitative effects of a distant reflector

We now attempt to investigate more quantitatively the effects of a distant reflector on the relation between the spectral index $`\mathrm{\Gamma }`$ and the iron line equivalent width and reflection amplitude $`R`$.

As noted in the previous section, the $`R`$-$`\mathrm{\Gamma }`$ relation depends on the variability mode. If one had a very accurate measurement of the $`R`$-$`\mathrm{\Gamma }`$ correlation we could derive the variability mode required in order to produce the observed correlation with a DCM. However due to the large uncertainties in the measurement of $`R`$ and $`\mathrm{\Gamma }`$ this is of little interest. We thus adopted the opposite approach: given the main observed characteristics of the variability in Seyfert 1s we test wether the presence of DCM can produce a correlation that is qualitatively similar to that observed. As discussed above, it seems that the spectral variability, in general, can be, at least roughly, described in term of a pivoting spectrum. We will thus assume that this pivoting spectrum represent the general (or average) variability mode of the Seyfert galaxies. We caution however that the details of the variability mode in Seyferts is not known, and the single pivoting point is only a simplifying approximation.

### 3.1 Model with constant cut-off energy

We assume that the time average primary spectrum seen by the DCM can be represented by a power law with a photon index $`\mathrm{\Gamma }=1.9`$ – corresponding to the average $`\mathrm{\Gamma }`$ in Seyferts 1 galaxies (Nandra et al. 1994; Matt 2001) – and an exponential cut-off at 400 keV. Neglecting any disc reflection, we further assume that the DCM geometry is such that the time averaged primary spectrum yield a reflection coefficient $`R=0.7`$ (close to the average observed $`R`$ for sources with $`\mathrm{\Gamma }1.9`$) . As a first order approximation, the shape of the reflection spectrum is computed using the pexrav procedure, i.e. assuming a slab reflector. Fixing the cut-off energy at 400 keV, we produced a set of instantaneous primary spectra for several photon indices $`\mathrm{\Gamma }`$ spanning the observed range 1.4–2.2. We did this for pivot points at energies 2, 5 and 10 keV respectively.

We then added the reflection component, produced as described above, to these instantaneous spectra, and fitted the resulting spectra with pexrav in the 2-30 keV range. We also fitted the simulated spectra in the 2–100 keV range. In both energy ranges, the derived best fit parameters were very similar. The $`R`$-$`\mathrm{\Gamma }`$ relations resulting from our modelling procedure are shown in Fig. 1 and compared with the correlation observed in the Ginga data.

Although the predicted correlation appears flatter than the observed one it *qualitatively matches the $`R`$-$`\mathrm{\Gamma }`$ correlation observed in the sample of Ginga data from ZLS99*.

### 3.2 Model with variable cut-off energy

Actually the amplitude of the reflection component is difficult to measure accurately and the derived $`R`$ values are quite dependent on the modelling of the primary emission (see Petrucci et al. 2001). In particular, the Ginga data used by ZLS99 to demonstrate the correlation do not enable one to constrain the cut-off energy. In their analysis they thus fixed $`E_c`$=400 keV in all their fits. Actually the $`R`$ derived from spectral fits depends on this assumed $`E_c`$ value. The BeppoSAX data allowed better constraints on the cut-off energy. Fitting these data with pexrav  shows that the cut-off energy differs from source to source and is also apparently correlated with the spectral slope. Fitting the $`E_c`$-$`\mathrm{\Gamma }`$ correlation presented by Perola et al. (2002), we found that this correlation can be qualitatively represented by the following function:

$$E_c=10^{0.73}\mathrm{\Gamma }^{5.7}\mathrm{keV}.$$
(1)
Perola et al. (2002) also suggested that the effects of the different cut-off energies could be responsible for part of the already mentioned differences between the BeppoSAX and Ginga$`R`$$`\mathrm{\Gamma }`$ correlations.

In order, to check the effect of a varying cut-off on the correlation produced by a DCM reflection, we generated a set of simulated spectra in a way similar to that described in section 3.1, but instead of a constant $`E_c`$, we introduced a “$`\mathrm{\Gamma }`$–dependent” $`E_c`$ as given by equation 1.

First, we mimicked the Ginga analysis of ZLS99 by fitting the set of simulated spectra in the 2-30 keV range assuming a constant cut-off at $`E_c`$=400 keV. As shown in the left panel of Fig. 2, this produces $`R`$-$`\mathrm{\Gamma }`$ correlations which are steeper than in the case of the models with constant $`E_c`$. In particular at low $`\mathrm{\Gamma }`$, the reflection amplitude is lower and reaches $`R`$=0 for the hardest spectra, as indeed observed by ZLS99. Thus allowing a correlation between $`E_c`$ and $`\mathrm{\Gamma }`$ results on a better agreement of the DCM model with the observations.

We then fitted the simulated spectra in the broader (2–200 keV) energy range that is allowed by BeppoSAX. The cut-off energy $`E_c`$ was then let as a free parameter and we consistently recovered best fit $`E_c`$ close to the injected ones. This procedure results in $`R`$-$`\mathrm{\Gamma }`$ correlations which are flatter and with higher $`R`$ values, consistent with the BeppoSAX correlation as can be seen from the right panel of Fig. 2.

We note that the $`E_c`$-$`\mathrm{\Gamma }`$ correlation found with pexrav in the BeppoSAX data does not necessarily reflect changes in the temperature of the Comptonising plasma. Indeed, although pexrav is useful to provide a simple and standard phenomenological description of observed spectra, this model represents only a rough approximation to real Comptonisation spectra. The physical interpretation of the fit parameters $`\mathrm{\Gamma }`$ and $`E_c`$ in terms of Thomson optical depth and temperature is not straightforward. Moreover the use of pexrav in spectral fits may produce or emphasize parameter correlations that are not necessarily found with realistic Comptonisation models (Petrucci et al 2000, 2001). To see if the $`\mathrm{\Gamma }`$-$`E_c`$ correlation obtained with pexrav could be produced at constant temperature, we used the compPS model111available at ftp://ftp.astro.su.se/pub/juri/XSPEC/COMPPS/ (Poutanen & Svensson 1996) to generate a set of realistic Compton spectra that we subsequently fit with pexrav. In the compPS simulations we assumed a spherical and isotropic geometry for the comptonising plasma with temperature fixed at $`kT_e`$=100 keV, and varied the Thomson optical depth of the comptonising plasma. We added to all spectra the same reflected flux. We fit these spectra with pexrav in the 2–200 keV range and indeed obtained a correlation between the resulting $`E_c`$ and $`\mathrm{\Gamma }`$ parameters. For $`\mathrm{\Gamma }<1.9`$, this correlation is however flatter than the observed $`E_c`$-$`\mathrm{\Gamma }`$ relation given by equation 1. The physical significance (if any) of the $`E_c`$-$`\mathrm{\Gamma }`$ correlation is thus unclear. From the phenomenological side however, this correlation provides the best pexrav representation of the data. It thus should be taken into acount when simulating observed spectra with this model.

We further note that the $`R`$-$`\mathrm{\Gamma }`$ correlation resulting from our compPS simulations with constant temperature, when these spectra are fit with pexrav in the Ginga range, is intermediate between what obtained for a constant $`E_c`$ (Fig. 1) and for $`E_c`$ given by equation 1 (left panel of Fig. 2).

### 3.3 The iron line equivalent width vs $`\mathrm{\Gamma }`$ relation

We now discuss the effects of a pivoting spectrum on the measured Fe line equivalent width ($`EW`$) produced by the DCM. For a power-law spectrum pivoting at an energy $`E_p`$, and a pure DCM spectrum, we have:

$$EW=\frac{I_{Fe}}{F_{E_p}}\left(\frac{E_p}{E_{K\alpha }}\right)^\mathrm{\Gamma },$$
(2)
with $`I_{Fe}`$ being the line photon intensity and $`F_{E_p}`$ the primary photon flux at energy $`E_p`$. From this relation it appears the $`EW`$ is positively correlated with $`\mathrm{\Gamma }`$ when the pivoting point is below $`E_{K\alpha }6.4`$ keV, while it is anti-correlated for higher pivoting energy (the $`EW`$ indeed represents the amount of reflection at 6.4 keV i.e. at a lower energy than the reflection bump). Of course, the normalization depends on the geometry of the DCM reflector. For an optically thick torus with a 30 $`\mathrm{deg}`$ opening angle, Ghisellini, Haardt and Matt (1994) derive an equivalent width $`EW80`$ eV. This estimate depends also widely on viewing angle and elemental abundances (Życki & Cerzny 1994; George & Fabian 1991). Using a large ASCA sample, LZ01 produced 3 average spectra of Seyfert 1s grouped according to increasing spectral index. For each spectra they could detect a narrow Fe line component that they attribute to reflection on a torus. Their estimates for the narrow line EW are shown in Fig. 3. It appears to be almost independent of $`\mathrm{\Gamma }`$ with an average EW of 50 eV.

Using a procedure similar to that described in section 3.1 and 3.2, we simulated the effects of spectral variability on the measured equivalent width (EW) of the iron line. We modeled the DCM iron line by a Gaussian with an intrinsic width $`\sigma =0.1`$ keV. The line flux was normalized to the observational result of LZ01, i.e. so that the $`EW`$ is 50 eV for the $`\mathrm{\Gamma }=1.9`$ and $`R=0.7`$ time-averaged primary spectra.

The resulting dependence of the measured equivalent width with $`\mathrm{\Gamma }`$ is plotted on Fig 3. We see that the trend (correlation or anti-correlation) in the $`EW`$-$`\mathrm{\Gamma }`$ relation is indeed qualitatively sensitive to the pivot point energy. However, if one consider that the pivot point location is likely to fluctuate (with time or from source to source), our assumed variability mode predicts essentially no specific relation between $`EW`$ and $`\mathrm{\Gamma }`$. Moreover, the fluctuations of $`EW`$ due to a pivoting spectrum are quite weak. Thus for a fixed geometry, our pivoting spectrum will produce an essentially constant $`EW`$, independently of $`\mathrm{\Gamma }`$, as observed by LZ01.

We note that the results of the averaging procedure of ZL01 should be taken with a caution (see Yaqoob et al. 2002). On the other hand, accurate measurements of the $`EW`$ of the narrow line in individual sources is somewhat difficult since it requires a high resolution observations and a good determination of the continuum and the broad-line component. Actually, the present estimates have large error bars and seem to depend on the modelling of the broad-line component and in particular, on the choice of the unconstrained parameters of the relativistic disc model. For instance, using XMM-Newton data of Mrk 359, O’Brien et al. (2001) find an $`EW`$ of the narrow line of 120 eV when assuming a disc extending down to 3 Rs, and $`EW=62`$ eV when fixing the disc inner radius at 50 Rs. Despite these uncertainties, the data on individual sources suggest that the $`EW`$ of the narrow line changes from source to source and/or on long time scale (i.e. longer than a few years). Although an anti-correlation between the narrow line $`EW`$ and source luminosity is reported (Nandra et al. 1997; Pound & Reeves 2002), the changes in the narrow line strength are apparently unrelated to $`\mathrm{\Gamma }`$.

Actually, the observed fluctuations of the narrow line $`EW`$ are likely to be dominated by intrinsic differences for instance in the geometry, elemental abundances or viewing angle of the reflecting material. Such differences are smeared out by the averaging procedure of LZ01 that shows the average $`EW`$-$`\mathrm{\Gamma }`$ (in)dependence.

Thus, our assumed pivoting spectrum induces only weak changes in the narrow line $`EW`$ that are observationally dominated by both uncertainties on the line strength measurement and a spread in the characteristics of the reflector from source to source. Simultaneously, the pivoting spectrum produces a strong dependence of $`R`$ on $`\mathrm{\Gamma }`$ that would dominate over the intrinsic differences from source to source and may contribute significantly to the reported $`R`$-$`\mathrm{\Gamma }`$ correlation. In this context, part of the observed spread in the $`R`$-$`\mathrm{\Gamma }`$ plane could be due to the slightly different DCM configurations in different sources. The prediction of the absence of strong correlation between the narrow line EW and $`\mathrm{\Gamma }`$ (linked with our assumed variability mode) is thus consistent with the present data and will be possibly tested when better data will be available.

## 4 The case of NGC 5548

From the previous section we see that the effects of DCM reflection, when important, make the geometrical interpretation of $`R`$ extremely misleading. *Trying to disentangle the temporal and geometrical effects is extremely difficult* and requires very long observations with time resolved spectral analysis. NGC 5548 was the target of such a monitoring campaign in 1999. It was observed in the X-rays on June 15-16, June 19-24, June 29-July 2, Jul 7 and August 16-18 with RXTE, thus providing a relatively good coverage of the source evolution during nearly 2 months. The results of these observations (as well as other contemporary multi-wavelength observations) are reported in Chiang et al. (2000, hereafter C00) . Interestingly, while the 2-10 keV flux changed by a factor of 2 during the campaign, C00 report a constant line flux suggestive of a DCM line. Strikingly the reflection amplitude $`R`$ was found to be essentially constant, with a trend for a correlation between $`R`$ and the 2-10 keV flux. The latter result suggests a close origin for the reflection bump and this is clearly inconsistent with DCM reflection. Then the line and the reflection bump should have a different origin, which appears difficult to explain.

In the case where the reflection component arises mainly from a distant reflector, the $`R`$ parameter is not very convenient to provide a physical insight. Rather the absolute reflection flux may be more relevant to probe the distant reflector since it is expected to be constant on short time scales (a few months) independently of the source flux and spectrum. On the other hand, in the case of close-by reflection, the reflected flux is expected to follow the fluctuations of the primary flux. We thus reanalyzed the C00 data, with a model of constant reflection spectrum using as fitting parameter, the absolute reflection flux $`RF`$ instead of the usual reflection amplitude $`R`$.

We used the 10 RXTE spectra taken at different epochs described in C00. Instead of using directly pexrav as did C00, our model consisted in an e-folded power-law (cutoffpl) plus a reflection component with a fixed shape (i.e. independent of the primary spectrum parameters). This reflection spectrum was produced by pexrav for a spectral index $`\mathrm{\Gamma }=1.9`$, a cut-off at 200 keV and standard abundances. The absolute normalization of the reflection component was a free parameter. As the reflection spectral shape is identical in all fits, this parameter is directly proportional to $`RF`$. The RXTE data do not enable one to constrain the high energy cut-off and we had to keep it frozen at an arbitrary value. We performed three fits of each spectrum assuming $`E_c`$=100, 200 and 400 keV respectively.

The dependence of $`RF`$ on the 2-10 keV flux is displayed on Fig. 4. For a fixed value of the cut-off energy the data are formally consistent with a constant value of $`RF`$. When the RF vs 2-10 keV flux relation is fit with a constant, the reduced $`\chi ^2`$ values are all below unity and the $`\chi ^2`$ probability is systematically larger than $`0.5`$. However, when a linear relation between the reflection flux and the 2-10 keV flux is assumed, we systematically get an improvement of the fit significant at better than the 95 per cent confidence level according to $`F`$ test. This correlation between the reflection normalization and 2–10 keV flux (in agreement with the finding of C00) is suggestive of a close reflector where the reflection flux increases almost linearly with the primary X-ray flux.

Nevertheless, as can be seen on Fig. 4, the measured $`RF`$ appears to depend on the cut-off energy. For our three fixed $`E_c`$ values we find 3 best fit $`RF`$ differing substantially. The $`RF`$ is systematically larger for lower $`E_c`$, the typical differences being about 40 per cent between $`E_c=400`$ keV and $`E_c=100`$ keV. Now, if we consider the uncertainties on the cut-off energy it appears that the observed $`RF`$ is also consistent within the error bars with a DCM reflection (i.e. constant $`RF`$). In particular one may consider that the actual $`E_c`$ might be positively correlated with the 2–10 keV flux.

Actually, the effects of the assumed $`E_c`$ on the $`RF`$ could be an artifact of our model associated with the particular e-folding cut-off shape assumed which may differ from a real thermal Comptonisation cut-off. We thus also model the primary emission using the more physical Comptonisation model compPS. We used the Thomson optical depth as the fitting parameter. Each spectrum was fit successively with the hot cloud temperature frozen at three different test values (40, 160, and 300 keV). On average, we obtain significantly lower $`RF`$s with the comptonisation model. The determination of the reflection is indeed very sensitive to the modelling of the primary continuum (see Petrucci et al. 2001). More precisely, as long as the temperature was fixed to a low value ($`<`$160 keV) the best fit $`RF`$ appears to be independent on the assumed temperature. There is a correlation between $`RF`$ and 2-10 keV flux (significant at the 98 per cent confidence level) consistent with a close reflector, although the DCM can not be formally ruled out. On the other hand at higher temperature it appears that $`RF`$ becomes strongly dependent on the assumed $`T_e`$. For $`T_e`$$`=`$300 keV we found significantly lower values for the $`RF`$ and almost no-correlation with the 2-10 keV flux. The improvement of the fit between a constant and linear model is significant only at the 63 per cent confidence level. This difference with respect to the low temperature case, is due to the fact that at large temperature the best fit optical depth is in general quite low ($`\tau _T0.2`$). The Comptonisation spectrum then differs from a pure power-law and this affects the measured $`RF`$.

The main conclusion here is that whatever the fitting model used, the measured reflection amplitude is very dependent on the assumptions on the cut-off energy (or temperature) which can be constrained, together with the reflection component, only with broad band instrument such as BeppoSAX or the future INTEGRAL.

As shown above, when these uncertainties are considered the present data become consistent with the DCM model. In particular, our results are fully consistent with a DCM reflector in the case of a high temperature Comptonising plasma, or a cut-off energy positively correlated with the 2-10 keV flux. It is worth noting that such high temperatures of order of 300 keV are indeed deduced from spectral fits to the the BeppoSAX data with realistic Comptonisation models (Petrucci et al. 2000) and that a correlation between the cut-off energy and spectral slope is inferred from the BeppoSAX sample using pexrav (see Matt 2001; Perola et al. 2002).

Independently, the narrow line component detected by Chandra (Yaqoob et al. 2001), recently confirmed by XMM-Newton (Pounds & Reeves 2002), together with the constant line flux reported by Chiang et al. (2000) are suggestive of DCM dominated reflection in this source. Thus, although the close-by reflection model is not ruled out by the Compton bump data alone, the DCM picture appears favored when these data are considered together with the iron line information.

## 5 Discussion

We considered the effects of reflection at large distance on the measured reflection amplitude $`R`$ and equivalent width of the Fe fluorescence line. We show that the remote cold material may in principle account simultaneously for the observed correlation between $`R`$ and $`\mathrm{\Gamma }`$, and the observed narrow line component which is predicted to be independent of $`\mathrm{\Gamma }`$.

We considered the $`R`$-$`\mathrm{\Gamma }`$ relation produced for different fixed pivoting point energies, with and without varying the cut-off energy, in all cases a correlation is produced. Actually, in the case of a constant comptonised flux and varying soft photon input, the pivot location is likely to fluctuate from the fact that changes in the cooling also affects the cut-off energy, and thus the energy range in which the comptonised power is released. As a consequence the pivot point may depend on $`\mathrm{\Gamma }`$. Note that if the spectral pivot shifts above 10-20 keV this produces an anti-correlation. Such an anti-correlation is reported by Larmer et al. (2000) in NGC 5506 where they find evidence for DCM reflection. However, even if the pivoting point energy changes among different sources or within the spectral evolution of an individual source we will just follow different correlations in the $`R`$-$`\mathrm{\Gamma }`$ plane and generally, a correlation will be produced. The correlation can be destroyed only if in a large fraction of the sources the pivot energy decreases very quickly with $`\mathrm{\Gamma }`$, or the pivot point is above a few 10 keV, in contrast to what is commonly observed.

Thus, despite the uncertainties a correlation between $`R`$ and $`\mathrm{\Gamma }`$ is expected in individual sources and is likely to be preserved in sample of sources. In the latter case this interpretation requires that most of the Seyfert sources present a similar variability mode (i.e. pivoting primary power-law), time average properties. as well as geometry of the DCM. While a pivoting powerlaw seems to be a general characteristics of Seyfert 1s, on the other hand, the geometry of the DCM is unknown. According to the unified model, the geometry of the DCM should be similar in all sources. In this context the small differences in geometry or variability mode would explain part of the spread of the data points in the $`R`$-$`\mathrm{\Gamma }`$ plane.

Our simulated $`R`$-$`\mathrm{\Gamma }`$ correlation are qualitatively similar to the observed one, but not exactly identical. In particular it seems that many of the data points lie below the predictions in both Fig. 1 and 2. This is not a fundamental problem for the DCM interpretation. This is simply a confirmation that the simple pivoting powerlaw model is a poor approximation to the actual variability mode. For instance if the pivoting points depends on Gamma and if for some reasons the pivoting point shift toward energies below 2 keV for the hardest spectra, it can produce as many sources with low reflection as required. However if we also consider the uncertainties on the measured $`R`$ this would be certainly overinterpreting the data. We thus prefered using a simple and approximative model that reproduces qualitatively the correlation, rather than a complicated and rather had hoc model.

In this study we neglected the contribution of reflection from the central part of the accretion flow. Actually the situation is probably more complex: the observed iron line profiles often present a broad component which is unlikely to be produced far away from the central engine. This broad component is probably associated with a inner reflection bump that may contribute to the $`R`$-$`\mathrm{\Gamma }`$ correlation. It is however not clear if the broad iron line component $`EW`$ is correlated with $`\mathrm{\Gamma }`$ as would be expected in such case (see LZ01 and Yaqoob et al. 2002 for a critical reevaluation of their results). At least, in several sources the broad line EW appears to be unrelated to the changes in $`\mathrm{\Gamma }`$ and present complicated behavior (Nandra et al. 1999, Vaughan & Edelson 2001 and references therein).

Another interesting point is the the analogy with X-ray binaries where a similar $`R`$-$`\mathrm{\Gamma }`$ relations is observed. Due to the time-scale involved this correlation cannot be accounted for by a DCM reflector (Gilfanov, Revnitsev & Churazov, 2000). The shape of the observed correlation seems however to differ between X-ray binaries and Seyferts, in particular the DCM can naturally explain both the very large and very low reflection amplitude that are not observed in X-ray binaries. Also unlike AGNs, several X-ray binaries present a *clear* iron line $`EW`$-$`\mathrm{\Gamma }`$ correlation. Actually accretion onto a black hole in an X-ray binary and an AGN is expected to proceed in a very different way due to the different boundaries conditions at large distance. The correlation in both types of sources could thus have a different origin.

In any case, the observations suggest that reflection arising from the DCM is important in many Seyfert galaxies. In particular, we showed that in NGC 5548 the reflection component could be dominated by the DCM reflector in the sense that, when uncertainties on the cut-off energy are taken into account, the RXTE data of C00 are consistent with a constant reflected flux. This would solve the previously reported inconsistency between the behaviour of the reflection bump and the observed constant line flux.

Finally we stress that our results do not exclude that part or even all of the $`R`$-$`\mathrm{\Gamma }`$ correlation in Seyferts is due to some contribution from the central region, however they show that this is not required. This has important implications for the modeling of the inner part of the accretion flow. It demonstrates that in Seyfert galaxies, the $`R`$-$`\mathrm{\Gamma }`$ correlation is not a good tool to discriminate between models for the physics and geometry of the X-ray emitting region. In particular emission models that do not produce a $`R`$-$`\mathrm{\Gamma }`$ correlation are not ruled out (see e.g. Malzac & Celotti 2002).

## 6 CONCLUSIONS

In this paper we attempted to reconcile the numerous evidences for DCM reflection in Seyfert 1s galaxies with the presence of a correlation between the reflection amplitude $`R`$ and the spectral index $`\mathrm{\Gamma }`$. We found that in the case of reflection on a DCM, the generally observed variability mode in Seyfert 1s naturally leads to a $`R`$-$`\mathrm{\Gamma }`$ correlation. We also showed that the measurement of $`R`$ and the shape of the observed $`R`$-$`\mathrm{\Gamma }`$ correlation is sensitive to the value of the cut-off energy. As long as the cut-off energy cannot be constrained, this effect significantly increases the error bars on the measured $`R`$ and possibly reduces the significativity of the correlation. In our case study of NGC 5548, we showed that when the effects of the cut-off energy are taken into account the data are fully consistent with the DCM picture, while in the context of close-by reflection the reported constant line flux remains a problem.

The EW of the narrow iron $`K_\alpha `$ line produced by reflection on the DCM is predicted to be generally independent of $`\mathrm{\Gamma }`$. This appears to be consistent with the present data, although the large uncertainties on the measurement of the the narrow line prevent us from reaching a firm conclusion.

In general, the possible presence of DCM reflection makes any attempt to constrain the physical properties and the geometry of the central engine, on the basis of the present reflection data, quite unreliable. Further observational studies, with instruments having a broad band coverage (INTEGRAL) and high resolution around the iron line (Chandra, XMM-Newton and latter XEUS), should allow one to disentangle the contribution of the DCM from that of the inner disc, and provide important insights on both the configuration of the DCM and the physics of the central engine.

## Acknowledgments

JM acknowledges grants from the Italian MURST (COFIN98-02-15-41) and the European Commission (contract number ERBFMRX-CT98-0195, TMR network ”Accretion onto black holes, compact stars and protostars”). We thanks James Chiang who provided us with the NGC 5548 data as well as Andrzej Zdziarski and Cesare Perola for the $`R`$-$`\mathrm{\Gamma }`$ data. We are also indebted to Laura Maraschi for useful comments on the manuscript.


text:
# Comment on the paper: ”Water content and its effect on ultrasound propagation in concrete- the possibility of NDE”. Ultrasonic 38(2000) 546-552 by Etsuzo Ohdaira and Nobuyoshi Masuzawa

## I Introduction

Concrete is one of the most important materials used is the development of a country’s infrastructure. Each year, thousands of tons of concrete is used in the restoration and rehabilitation of buildings. The development of the construction industry has helped in finding more effective and economic methods of inspection, which permit quicker and more reliable control. For this reason, the application of non destructive test methods is becoming more widely used in the industry, among then the ultrasound Ref. Nevertheless, the effective application of the method in the study of concrete depends on the knowledge of a series of factors that affect the pulse velocity as well as the mechanical parameters Ref. One of these factors is the water content of the concrete , in terms of which the ultrasonic pulse velocity can increase or decrease significantly , while the variations in the compressive strength remains insignifican Ref. The incidence of this factor can bring about imprecision in the estimation of the mechanical properties , specifically the strength of the concrete. The results obtained in Refare being analyzed in this paper and a model adequately describing the behavior of the velocity in terms of the water content is presented.

## II Methods of measurement

In the work Ref, cylindrical samples of 200x100mm (diameter)are used for various dosifications. The specimen is situated on the scale to measure the variation of the water content according to weight lost ( until weight is constant). The transducer are attached using the transmission method, which allows the measurement of the velocity in terms- of the water content.

## III Discussion of the results.

In Ref, it is proposed that behaviour of the velocity as function of the concentration of water shows a lineal behavior, which is not particularly precise. If the linnear model is considered , it must be admitted that an increase in the concentration, in the absorption process, would always provoke an increase in velocity which would not always be possible from the physical point of view because the saturation process is asymptotic, meaning that the velocity must reach an asymptotic maximum. In figures 1,2,3, the non linnear tendency of the process is shown. From the calculation of the probability of the goodness of fit it is observed that in all cases , there are no significant differences amoung velocity for a significance level of 5

It is then necessary to find a model which adequately describes the phenomenon under study. In spite of the fact that a 3 rd order polynomial describes this process rather well from a statistical point of view, it lacks common sense when seen from a physical point of view. We have a phenomenon in which the extreme values have an asymptotic behavior as a result of the drying process of the sample. As it is known, the increase of water in the material results in the increase in the velocity of ultrasonic signal. So, if it is taken into account that the material has a saturation point which is reached asymptotically, then the velocity also has a maximum value reached asymptotically. The former considerations can be written as follows:

Although the decrease in the concentration is proportional to the decrease in velocity, this decrease is greater while the velocity moves further away from a certain maximum value

This affirmation can be expressed by the following differential equation:

$$\frac{dV}{dC}=k(V_0V)$$
(1)
And the solution is:

$$V=V_0A_0e^{kx}$$
(2)
Where: the integration constant A0 is considred negative in terms of the observations made the experiment. v-velocity V0, A0 parameters of the model C-water content K-parameter of the model

Adjusting for each one the proportions , the following coefficients are obtained

| Parameters | Test Piece A-1 | Test Piece B-1 | Test Piece C-1 |
| --- | --- | --- | --- |
| $`V_0`$ | $`4487.58235`$ | $`4406.20072`$ | $`4234.08337`$ |
| $`A_0`$ | $`553.01462`$ | $`523.59085`$ | $`429.80387`$ |
| $`\frac{1}{k}`$ | $`2.88231`$ | $`2.20381`$ | $`3.07733`$ |
| $`R^2`$ | $`0.97`$ | $`0.97`$ | $`0.92`$ |

In all cases , the models converge rapidly.

From (2) it is observed that as concentration decreases, velocity decreases gradually and as C tends to zero , velocity assumes its minimum value, analogously , as C tends to infinite, velocity achieves its maximum value, describing a gradual process which corresponds with the capillary desorbption and diffusion phenomenon, which depends on time and concentration gradient. There are two stages involved in this drying process: initially capillarity, where the free water molecules rise to the pores surface and evaporate into the surrounding air, and the diffusion period, which is responsible for the evaporation of water trough the fine pores. In order to know how capillarity varies with time and concentration, it is necessary to solve the following equations for specific boundary condition of the desorbption process Ref

$$\frac{\theta }{t}=D\theta $$
(3)
Where D-capillary diffusivity t-time $`\theta `$— concentration of water.

Capillarity as well as diffusion, depend in great measure on the porosity of the medium. Therefore the coefficients are specific for each type of concrete.

Although the loss of water accompanied by contractions corresponds to a cause and effect situation, their relationship is not a simple one. When concrete begins to dry up, the water lost first is the free water molecules contained in the capillaries , which practically produces no contraction Ref. For this reason , the first experimental point where a sudden variation of velocity occurs can be rejected.

Bearing in mind what has been said so far, an the existing dependence between velocity and the acoustic characteristic of the medium , it can then be understand why velocity varies whit the sample water content Ref. There is a significant difference between the acoustic impedance of the air which exists initially in the superficial pores and that of the water which is absorbed in the course of time , in other words , if the pores are filled with air , the medium is denser , which favors the propagation’s of the high- frequency components which are more sensitive to these change. As observed in Ref, on obtaining the frequency spectrum, the value of the fundamental component decreases(attenuate) as the concentration of water decreases. However, the acoustic impedance of water and air are small in comparison with impedance of concrete. Why is it then that this variations in velocity occurs? It is logical that the wave is propagated preferably through the solid matrix, but if the number of pores is very high, then velocity is affected by the presence of pores. Besides, water penetrates troungh a distance of 3-4 cm in the interior of the sample. This implies that all the area near the surface is saturated and its dimensions are of the order of the wavelength of the fundamental component , meaning that for velocities above 4000 m/s and frequencies of the order of 55- 100 KHz , the wavelength of the fundamental component is about 4-6 cm, the entire pulse front being affected by the presence of water.

It is important to demonstrate that the relative error of estimated velocity can be calculated from the model for any given concentration using the following expression:

$$\delta V=0.0108895+\frac{0.966460.0092429C}{0.123232+e^{0.3469446C}}$$
(4)
where: delta V——– relative error of velocity C ——concentration.

This error, as observed in Figura 4, it is approached to one percent in the low water content zone and diminishes with the increase of concentration as a consequence of the asymptotic tendency of the saturation process.

It is known that the fundamental objective of the application of ultrasound in the study of concrete is to correlate velocity witch compressive strength Ref, however when the humidity condition changes as has been analyzed, the behavior of the strength is not analogous to that of velocity, it practically maintains an almost constant value, as shown in Ref for this reason, it is very important to keep in mind that high velocity do not always imply that values of strength will also be high, so mistakes can be made in the estimation of strength. In Figure 5 it is shown the tendency of estimated strength using the expression (5) Ref, employing the value of experimental velocity obtained in Ref, the velocity estimated using the linnear model proposed in Ref and velocity estimated through the exponential model proposed in this work. In all cases velocity increases with water content. It is observed that estimated strength also augments with water content, which is incorrect. That’s why to estimate the strength from velocity it is necessary to know moisture percent of the sample.

$$R=aVb$$
(5)
where: R- compressive strength (Mpa) a,b-parameters of the model V-velocity

In this case we employ a=0.0447 y b= 149.67 Ref

## IV conclusions

It is proven physically and statistically that the model which best describes the behaviour of velocity with respect to water concentration is the decreasing exponential model, from which the compressive strength can be estimated with a good confidence.